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SUMMARY 

The  main  objective  of  this  research  is  the  development  and  verification  of  numeri¬ 
cal  models  for  simulating  the  dynamic  and  thermodynamic  characteristics  of  the 
Arctic  ice  cover.  Emphasis  has  been  placed  on  realistically  parameterizing  the  ridg¬ 
ing  process  within  the  framework  of  a  variable  thickness  sea  ice  model. 

The  basic  model  consists  of  a  multiple-level  ice  thickness  distribution  coupled  to  a 
complete  momentum  balance  employing  a  viscous  plastic  constitutive  law.  A  flow 
chart  describing  the  overall  model  is  given  in  Figure  1.  Note  that  a  full  heat  budget 
code,  a  fixed  depth  oceanic  mixed  layer,  and  a  sea  ice  thermodynamic  model  are  used 
in  conjunction  with  the  ice  thickness  evolution  equations.  Details  of  the  equations 
are  presented  in  the  main  text.  In  the  numerical  code,  the  ridge  redistribution  pro¬ 
cess  is  formulated  in  a  general  manner  so  that  a  variety  of  different  ridging  mecha¬ 
nisms  may  be  modeled.  The  numerical  scheme  is  formulated  in  a  fixed  Eulerian  grid 
so  that  integrations  over  unlimited  time  intervals  may  be  performed.  In  addition,  in 
the  mixed  layer  formulation,  lateral  melting  terms  are  included.  In  the  computer 
code,  special  attention  has  been  given  to  conservation  properties.  Specifically,  the 
code  is  formulated  to  manifestly  conserve  air-sea  heat  exchanges  so  that  all  energy 
exchanges  are  precisely  accounted  for  in  terms  of  ice  formation,  ice  melt  or  warming 
of  the  oceanic  mixed  layer.  Because  of  these  conservation  properties,  the  code  is 
suitable  for  coupling  to  atmospheric  and  oceanic  general  circulation  models. 

Because  of  the  importance  of  ridging  in  this  model,  considerable  effort  has  gone 
into  designing  a  realistic  ridge  distribution  process.  To  this  end  a  redistribution  pro¬ 
cess  consistent  with  observed  and  hypothesized  physics  of  the  ridging  process  is  pro¬ 
posed.  This  redistributor  is  analytically  examined  and  compared  to  a  previous  pro¬ 
posed  redistributor  (Fig.  3  and  4).  The  two  key  physical  features  incorporated  in  this 
redistributor  are  1)  under  deformation,  ridging  transfers  ice  to  a  variety  of  thickness 
categories,  2)  typical  ridge  heights  increase  more  slowly  than  does  the  thickness  of 
ice  being  ridged.  Analytic  calculations  show  that  taking  into  account  the  second  fea¬ 
ture  tends  to  preferentially  increase  thin  ice  strengths  as  compared  to  previously  pro¬ 
posed  redistributors  (Fig.  4). 

The  main  approach  to  analyzing  the  fully  coupled  model  characteristics  has  been 
to  carry  out  an  equilibrium  variable  thickness  simulation  of  the  Arctic  Basin  ice  cov¬ 
er.  In  this  simulation,  a  plastic  rheology  coupled  to  a  10-level  ice  thickness  distribu¬ 
tion  was  found  to  yield  realistic  geographical  variations  in  ice  thickness,  with  ice  in 
excess  of  7  m  thick  along  the  Canadian  Archipelago  and  about  2  m  along  the  Alas¬ 
kan  North  Slope  (Fig.  7).  Examination  of  this  buildup  shows  that  it  takes  several 
years  to  fully  evolve  and  is  largely  due  to  ridging. 

The  volume  of  ridged  ice  formed  each  year  also  shows  interesting  spatial  varia¬ 
tions  (Fig.  11).  The  general  shapes  of  the  roughness  contours  in  this  figure  agree  well 
with  surface  roughness  observations.  The  major  feature  is  a  heavy  buildup  of  ridg¬ 
ing  off  the  Canadian  Archipelago  together  with  less  ridging  near  the  pole.  In  addi¬ 
tion,  there  is  a  zone  of  heavy  ridging  just  off  the  North  Slope  that  is  in  agreement 
with  observations  there.  A  final  feature  is  a  tongue  of  high  ridging  further  offshore 
near  the  tip  of  Greenland.  Roughness  data  obtained  from  submarine  profiles  also 
show  such  a  tongue. 


Several  sensitivity  experiments  were  performed  to  determine  the  major  factors 
dictating  the  ice  edge  location  (Fig.  7,  13,  and  14).  The  results  indicate  that  it  is  a 
combination  of  summer  melt  coupled  with  large  amounts  of  offshore  advection  ear¬ 
ly  in  the  spring.  The  most  critical  factor,  however,  appears  to  be  summer  melt.  An 
increase  in  the  albedo  of  melting  ice  by  only  10%,  for  example,  substantially  im¬ 
proved  the  results  from  North  Slope  ice  edge  (Fig.  13). 

The  results  also  give  some  insight  into  strengths  needed  for  proper  simulations. 
Basically,  by  use  of  the  ridging  parameterization  mentioned  above,  a  realistic  ice 
buildup  was  simulated.  However,  the  simulated  ice  velocities  were  rather  large 
(Table  1).  Sensitivity  analyses  suggest  that  a  modest  increase  in  frictional  losses 
would  be  adequate  for  obtaining  more  realistic  ice  velocities  and  strengths.  Further, 
more  detailed  studies  with  polar  drifting  buoys  are  needed  to  more  precisely  deter¬ 
mine  these  strengths. 

With  regard  to  future  research,  substantial  progress  has  been  made  on  the  devel¬ 
opment  and  numerical  examination  of  a  complete  variable  thickness  sea  ice  model. 
However,  considerable  theoretical  and  numerical  work  remains  to  determine  more 
precisely  the  strengths  and  weaknesses  of  this  model  as  it  stands.  The  ice  dynamics 
formulation,  particularly,  needs  further  study. 
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NUMERICAL  MODELING  OF  SEA  ICE 
DYNAMICS  AND  ICE  THICKNESS 
CHARACTERISTICS 
A  Final  Report 

William  D.  Hibler  III 


INTRODUCTION 

In  the  polar  regions  the  interaction  between  the  atmosphere  and  the  ocean  is  significantly 
modified  by  the  presence  of  a  sea  ice  cover.  Typically,  this  ice  cover  contains  a  variety  of  ice 
thicknesses  that  evolve  in  response  to  both  dynamic  and  thermodynamic  forcing.  This  vari¬ 
able  thickness  feature  can  substantially  change  the  way  the  ice  cover  modifies  the  air-sea  in¬ 
teraction  and,  hence,  affect  the  atmospheric  and  oceanic  circulation.  Numerical  and  empiri¬ 
cal  studies  demonstrating  the  relevance  of  this  interaction  to  the  atmosphere  have  recently 
been  carried  out  by  Herman  and  Johnson  (1978),  Walsh  and  Johnson  (1979a),  Saltzman  and 
Moritz  (1980),  Lemke  et  al.  (1980)  and  Manabe  and  Stouffer  (1980).  In  light  of  these  results, 
properly  simulating  a  variable  thickness  sea  ice  cover  has  become  particularly  relevant  to  nu¬ 
merical  investigations  of  climate. 

Recent  research  on  modeling  a  variable  thickness  ice  cover  is  largely  based  on  pioneering 
work  by  Coon  (1974),  Thorndike  et  al.  (1975)  and  Rothrock  (1975).  Specifically,  Thorndike 
et  al.  (1975)  introduced  an  areal  ice  thickness  distribution  function  and  developed  equations 
for  the  dynamic-thermodynamic  evolution  of  this  distribution.  The  input  fields  for  the 
evolution  equations  consist  of  a  two-dimensional  ice  velocity  field  and  ice  growth  rates  ver¬ 
sus  thickness  and  time.  In  this  model  the  effects  of  ice  dynamics  on  the  ice  thickness  distribu¬ 
tion  are  quantified  by  allowing  thin  ice  formation  and  redistribution  of  ice  to  thicker  catego¬ 
ries  (because  of  ridging)  in  response  to  deformation.  Similarly,  thermodynamic  effects  cause 
a  rearrangement  of  the  relative  amounts  of  ice  in  different  categories.  Rothrock’s  (1975)  con¬ 
tribution  was  to  provide  a  means  for  coupling  the  ice  thickness  distribution  to  the  rheological 
behavior  of  sea  ice.  This  rheological  behavior,  in  turn,  affects  the  ice  dynamics.  For  this  pur¬ 
pose  Rothrock  (1975)  noted  that  the  rate  of  work  done  on  the  ice  through  ridging  is  related  to 
the  work  done  by  the  ice  interaction  forces.  By  combining  this  idea  with  the  concept  of  a 
plastic  constitutive  law  for  sea  ice  developed  by  Coon  (1974),  it  is  possible  to  form  a  coupled 
set  of  equations  describing  the  dynamic-thermodynamic  behavior  of  sea  ice. 


Numerical  investigations  of  the  variable  thickness  concept  have  mainly  concentrated  on 
examining  the  ice  thickness  distribution  independent  of  its  coupling  effects  on  ice  dynamics. 
Thorndike  et  al.  (1975),  for  example,  used  the  deformation  field  for  a  lagrangian  parcel  of 
ice  (defined  by  three  contemporaneous  drifting  stations  in  the  Arctic  Basin)  together  with 
growth  rate  estimates  versus  thickness  and  time  to  drive  a  thickness  distribution  model.  Sim¬ 
ulations  with  these  data  reproduced  an  ice  thickness  distribution  in  qualitative  agreement 
with  observations.  Quantitative  discrepancies  in  the  results  did,  however,  occur  in  the  thick 
end  of  the  thickness  distribution.  Subsequently,  Maykut  (1978)  used  observed  and  simulated 
thin  ice  percentages  to  estimate  the  regionally  averaged  heat  input  to  the  atmosphere.  His  es¬ 
timates  suggest  that  in  the  central  Arctic  in  winter,  heat  exchange  through  ice  in  the  0-1 .0  m 
thickness  range  is  approximately  equal  to  heat  exchange  through  thicker  ice.  This  is  so  even 
though  the  thicker  ice  constitutes  a  much  larger  areal  fraction  of  the  ice  cover. 

The  simulations  by  Thorndike  et  al.  (1975)  and  Maykut  (1978)  established  the  validity  and 
importance  of  the  ice  thickness  distribution  concept.  However,  they  were  limited  in  the  sense 
that  changes  in  the  ice  thickness  distribution  were  not  allowed  to  modify  the  ice  dynamics.  In 
addition,  a  lagrangian  element  that  covered  only  a  small  portion  of  the  basin  was  used.  Be¬ 
cause  of  this  localized  lagrangian  formulation,  relative  geographical  variations  were  not  ex¬ 
amined  and  advection  effects  were  not  explicitly  modeled.  A  useful  way  to  examine  the  be¬ 
havior  of  a  variable  thickness  dynamic-thermodynamic  sea  ice  model  that  is  more  fully  coup¬ 
led  is  to  carry  out  a  seasonal  equilibrium  simulation.  Specifically,  by  integrating  such  a  mod¬ 
el  over  sufficiently  large  time  intervals  (several  years),  results  for  both  drift  and  thickness  can 
be  obtained  that  are  relatively  independent  of  initial  conditions.  While  a  number  of  short¬ 
term  integrations  of  such  models  have  been  done  for  localized  regions  (Coon  et  al.  1976, 
Pritchard  et  al.  1977),  no  equilibrium  simulations  have  been  performed.  An  impediment  to 
conducting  such  a  fully  coupled  dynamic-thermodynamic  equilibrium  simulation  is  the  com¬ 
putational  difficulty  in  solving  the  equations.  In  particular,  carrying  out  coupled  numerical 
simulations  over  times  long  enough  to  examine  the  thermodynamic  evolution  of  the  ice  cover 
requires  explicit  inclusion  of  nonlinear  advection  terms.  Such  terms  are  not  included  in  either 
the  Thorndike  et  al.  (1975)  simulations  or  the  Coon  et  al.  (1976)  and  Pritchard  et  al.  (1977) 
studies. 

This  report  describes  a  numerical  framework  suitable  for  simulating  a  variable  thickness 
sea  ice  cover  over  a  seasonal  cycle  and  presents  a  new  redistributor  for  modeling  ridge  build¬ 
up  on  the  geophysical  scale.  In  addition  to  seasonal  simulations,  this  model  also  can  be  ap¬ 
plied  to  shorter-term  forecasts.  While  it  has  certain  features  in  common  with  the  model  de¬ 
veloped  by  Hibler  (1979)  in  an  earlier  study  of  the  Arctic  ice  cover,  the  present  framework 
contains  a  more  general  treatment  of  the  ice  strength,  the  ice  thickness  distribution,  and  the 
ice  growth  and  decay.  By  combining  this  framework  with  a  thermodynamic  sea  ice  model 
similar  to  that  of  Semtner  (1976),  a  seasonal  equilibrium  simulation  of  the  Arctic  Basin  ice 
cover  is  performed.  A  flow  chart  describing  the  overall  model  is  given  in  Figure  1. 

Several  shorter  sensitivity  studies  are  also  carried  out.  In  analyzing  the  results  of  these  sim¬ 
ulations,  particular  attention  is  paid  to  the  evolution  and  sensitivity  characteristics  of  the 
summer  ice  edge,  and  the  ice  thickness  buildup  through  ridging  off  the  Canadian  Archi¬ 
pelago.  These  characteristics  have  not  been  examined  in  previous  variable  thickness  studies. 
To  a  large  degree  this  study  was  motivated  by  the  lack  of  a  variable  thickness  equilibrium 
simulation  of  the  Arctic  Basin  ice  cover.  The  main  purpose  of  this  report  is  to  present  a  com¬ 
plete  variable  thickness  sea  ice  model  and  to  investigate  the  degree  to  which  such  a  model 
with  full  dynamic-thermodynamic  coupling  can  reproduce  the  observed  thickness,  drift  and 
ridging  characteristics  of  the  Arctic  Basin. 
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Figure  1.  Flow  chart  for  variable  thickness  dynamic-thermodynamic 
sea  ice  model.  C— Coriolis  force;  tw—  water  stress  due  to  ice  motion; 
ta— air  stress;  F— internal  stress  variation;  G — ocean  currents;!— ocean  tilt; 
u — ice  velocity;  Djj/Dt — ice  acceleration  and  momentum  advection;  m — ice 
mass  per  unit  area;  g—ice  thickness  distribution;  f(h) — ice  growth  rate; 
h— ice  thickness;  \pfg,h)—ice  thickness  redistribution. 


DESCRIPTION  OF  THE  MODEL 

A  fully  coupled  dynamic-thermodynamic  sea  ice  model  can  be  divided  into  the  following 
components:  a  momentum  balance  describing  ice  drift,  which  includes  air  and  water  stresses, 
Coriolis  force,  internal  ice  stress,  inertial  forces  and  ocean  tilt;  an  ice  rheology,  which  relates 
the  ice  stress  to  the  ice  deformation  and  strength:  ice  thickness  distribution  equations,  which 
describe  the  evolution  of  the  ice  thickness  characteristics  caused  by  thermodynamic  and  dy¬ 
namic  effects;  and  an  ice  strength,  determined  as  a  function  of  the  ice  thickness  distribution. 
An  additional  component  needed  for  long-term  integration  of  the  model  is  a  thermodynamic 
code,  which  specifies  the  growth  and  decay  rates  of  various  ice  thicknesses  in  the  environ¬ 
ment  of  a  sea  ice  cover  of  variable  thickness.  The  first  two  components  are  discussed  by 
Hibler  (1979)  in  a  paper  describing  a  dynamic-thermodynamic  sea  ice  model  with  two  thick¬ 
ness  levels.  In  this  model  the  dynamical  equations  include  air  and  water  stresses,  Coriolis 
force,  internal  ice  stress,  inertial  forces,  momentum  advection  terms  and  ocean  tilt.  Non¬ 
linear  boundary  layers  for  both  the  ocean-ice  and  air-ice  surface  traction  are  used.  For  the 
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Figure  7  (cont  ’d). 


reaches  its  maximum  and  minimum  extremes  in  thickness  and  extent.  For  comparison  the 
August  plots  show  observed  ice  edge  estimates  obtained  from  the  British  Meteorological  Of¬ 
fice  at  Bracknell.  This  ice  edge  is  the  mean  40%  ice  concentration  at  the  end  of  August  over 
the  time  period  1966-1976.* 

The  dominant  contrast  between  the  two  cases  is  a  pronounced  difference  between  the  geo¬ 
graphical  variations  of  the  ice  thicknesses.  The  results  considering  only  thermodynamic  ef- 


1  Personal  communication  with  G.  Rowland,  British  Meteorological  Office,  Bracknell,  England,  March  1980. 


A  125-km  resolution  grid  was  used  for  the  simulations  and  is  shown  in  Figure  5  (this  grid  is 
identical  to  that  used  by  Hibler  [19791).  Some  of  the  numbered  grid  cells  are  referred  to  later 
in  the  text.  The  shaded  grid  cells  near  Greenland  are  taken  to  be  “outflow”  grid  cells.  Ice  is 
only  allowed  to  be  transferred  into  these  grid  cells  by  advection,  and  once  there  is  considered 
to  have  “flowed  out”  of  the  basin.  More  details  on  the  treatment  of  ice  dynamics  and  advec¬ 
tion  at  such  cells  are  given  in  Hibler  (1979). 

For  the  variable  thickness  simulations,  the  key  results  of  the  heat  budget  calculation  are 
the  ice  growth  rates  versus  thickness.  Figure  6  shows  the  seasonal  growth  rates  (plotted  every 
16  days)  for  different  thicknesses  at  grid  cell  8.  Also  plotted  is  the  wind  speed  at  the  grid  cell. 
The  main  feature  is  a  pronounced  seasonal  cycle  in  growth  rates  that  becomes  more  accentu¬ 
ated  for  thinner  ice  types.  Superimposed  on  the  seasonal  cycle  are  fluctuations  in  growth 
rates  attributable  to  varations  in  the  sensible  and  latent  heat  losses  which  depend  heavily  on 
the  wind  speed.  This  heavy  dependence  of  the  thin  ice  and  open  water  growth  on  wind  speed 
has  been  experimentally  observed  by  Andreas  (1980)  and  is  a  dominant  feature  of  Maykut’s 
(1978)  numerical  heat  exchange  results  for  thin  ice.  Note  that  under  summer  conditions  the 
melting  rate  of  ice  is  not  dependent  on  the  thickness  because  of  use  of  a  fixed  albedo  for 
melting  ice.  However,  while  the  total  melting  rate  is  the  same,  the  thinner  ice  will  have  a 
greater  ratio  of  bottom  melt  to  top  melt  owing  to  the  more  efficient  conduction  of  heat  into 
the  mixed  layer  (see  Appendix  B). 

Basin-wide  ice  thickness  and  velocity  characteristics 

After  5  years  both  simulations  approach  a  seasonal  equilibrium,  with  the  ice  thickness  and 
(in  the  standard  case)  ice  velocity  characteristics  changing  little  between  corresponding  days 
of  sequential  years.  The  basin-wide  thickness  characteristics  for  the  fifth  year  of  integration 
are  shown  in  Figure  7  for  April  and  August.  During  these  months  the  pack  ice  approximately 
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Figure  7.  Average  April  and  August  thickness  contours  (m) 
for  the  fifth  year  of  the  standard  and  "thermodynamics- 
only”  simulations.  The  dashed  line  represents  the  average  40 % 
:ce  concentration  contour  at  the  end  of  August  for  the  time  period 
1966-1976  as  estimated  by  the  British  Meteorological  Office, 
Bracknell,  England. 
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creased  to  0.66  (from  0.616  in  the  standard  case)  and  the  oceanic  heat  flux  was  set  equal  to 
zero.  In  the  “strength”  sensitivity  study  the  frictional  losses  in  ridging  were  taken  to  be  nine 
times  the  potential  energy  change  (i.e.,  C  =  10  C„  in  eq  3  of  Appendix  A). 

Input  fields  to  drive  the  model  consisted  of  monthly  climatological  air  temperatures  and 
dew  points  (from  Crutcher  and  Meserve  1970),  together  with  observed  winds  that  were  aver¬ 
aged  over  8  days  and  varied  with  time  over  the  year-long  period:  May  1962  to  May  1963. 
(This  particular  time  was  chosen  because  of  the  simultaneous  presence  of  one  U.S.S.R.  and 
two  U.S.  drifting  ice  stations  that  provided  observed  ice  drift  data).  For  the  calculation  of 
geostrophic  ocean  current  fields,  mean  dynamic  topography  values  reported  by  Coachman 
and  Aagaard  (1974)  were  used.  For  the  radiation  code,  parameterizations  similar  to  those 
employed  by  Parkinson  and  Washington  (1979)  were  used.  Specifically,  daily  global  solar  ra¬ 
diation  under  cloudless  skies  was  obtained  by  integrating  an  empirical  equation  by  Zillman 
(1972)  over  solar  zenith  angles  for  any  particular  day.  (Zenith  angles  at  half-hour  intervals 
for  this  purpose  were  obtained  from  a  numerical  solution  of  Kepler’s  equation.)*  Incoming 
longwave  radiation  was  obtained  using  Idso  and  Jackson’s  (1969)  formula  for  clear  skies. 
For  cloud  cover  estimates,  values  from  Huschke  (1969)  as  reported  by  Parkinson  and  Wash¬ 
ington  (1979)  were  employed.  Readers  interested  in  more  details  on  the  various  climatologi¬ 
cal  and  radiation  forcing  fields  are  referred  to  Parkinson  and  Washington  (1979). 


Figure  6.  Growth  rates  versus  time  for 
different  ice  thickness  categories  at 
grid  cell  8.  Also  plotted  for  comparison  is 
the  wind  speed  versus  time  at  this  location. 
Data  points  are  plotted  every  16  days  and 
the  ice  thicknesses  are  labeled  in  metres. 


'Personal  communication  with  H.  Holloway,  Geophysical  Fluid  Dynamics  Laboratory,  Princeton,  New  Jersey, 
November  1979. 


Figure  5.  Fixed  square-mesh  grid  used  for  numerical  calculations.  The  num¬ 
bered  grid  cells  denote  locations  where  time  series  of  ice  characteristics  are  moni¬ 
tored. 


To  investigate  the  behavior  of  this  fully  coupled  variable  thickness  sea  ice  model,  several 
seasonal  simulations  of  the  Arctic  Basin  ice  cover  were  carried  out.  These  simulations  were 
performed  by  combining  the  variable  thickness  equations  developed  here,  i.e.,  eq  4-7  and  eq 
A5,  A7,  A8  and  A10,  with  a  previously  developed  viscous-plastic  ice  dynamics  model  (eq  1- 
1 1  of  Hibler  [1979]).  For  the  standard  experiment  the  resulting  dynamic-thermodynamic  sys¬ 
tem  of  equations  was  numerically  integrated  for  5  years  at  1-day  time  steps  using  forcing 
fields  with  a  1-year  periodicity.  This  integration  time  was  found  to  be  adequate  to  obtain  sea¬ 
sonally  varying  equilibrium  results.  For  comparison,  the  thermodynamic  equations  only 
were  also  integrated  for  5  years  with  the  same  forcing.  In  both  cases  the  integration  was 
started  on  1  January  with  mean  ice  thicknesses  of  3.2967  m  ( =  3.0  x  10!  kg  m'Vg,  where  q,  is 
the  density  of  ice)  at  all  grid  cells. 

In  the  standard  case  this  mean  thickness  was  produced  by  including  ice  in  the  three  thick¬ 
ness  categories  centered  at  1 .46,  2.61  and  4.23  m  (see,  e.g.,  Fig.  2).  The  wind  and  water  drag 
coefficients  and  Coriolis  parameter  were  identical  to  those  used  by  Hibler  (1979).  The  re¬ 
maining  constants  are  the  mixed  layer  depth,  dmi ,  (set  at  30  m  following  Semtner  [1976]),  and 
various  thermodynamic  parameters  (viz  the  surface  albedo  of  ice  under  melting  and  freezing 
conditions,  the  open  water  albedo,  the  oceanic  heat  flux,  and  sensible  and  latent  heat  trans¬ 
fer  coefficients).  The  values  for  these  thermodynamic  parameters  for  the  standard  case  are 
given  in  Appendix  B.  In  addition  to  the  standard  case,  two  shorter  1-year  dynamic-thermo- 
dynamic  simulations  were  carried  out  using  the  fourth  year,  31  December,  equilibrium 
simulation  results  as  initial  conditions.  The  purpose  of  these  shorter  simulations  was  to 
assess  the  sensitivity  of  the  equilibrium  results  to  ice  strength  and  thermodynamic  param¬ 
eters.  In  the  “growth”  sensitivity  simulation,  the  ice  albedo  under  melting  conditions  was  in- 
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Figure  4.  Ice  strength  versus  thickness  for 
different  redistributors.  The  solid  and  long- 
dashed  lines  are  for  a  uniform  redistributor 
with  the  maximum  thickness  scaling  as  the 
square  root  (Hibler)  of  the  thickness  of  ice  be¬ 
ing  ridged  (K  =  Vfu),  and  linearly  (modified 
Hibler)  with  the  thickness  of  ice  being  ridged 
( K  =  15).  Frictional  losses  are  assumed  to  be 
equal  to  changes  in  gravitational  potential 
energy. 


P *  = 


Chi  (4P-3P-1) 

3(r-  or 


where  f  =  -J H*/h0  for  the  Hibler  redistributor  and  f  =  K/ 2  for  the  modified  Hibler  case, 


P*  =  ChlK 


(Thorndike  et  al.). 


P*  =  Ch,H’ 


(rubble  redistribution).  With  a  representative  value  for  C  of  0.8026  x  10"’  N  m"\  the  various 
scalings  in  Figure  3  yield  the  strength  versus  thickness  results  shown  in  Figure  4.  In  the  con¬ 
stant  maximum  cutoff  case  we  have  used  the  modified  Hibler  redistribution.  However,  from 
the  above  strength  equation  it  is  clear  that  the  Thorndike  et  al.  definition  would  give  about 
the  same  strength  versus  thickness  variation,  but  with  al)  strengths  50%  larger  than  the  modi¬ 
fied  Hibler  redistribution. 

The  main  feature  illustrated  by  Figure  4  is  that,  as  expected,  the  square  root  scaling  tends 
to  give  higher  strengths  for  thinner  ice  than  the  constant  scaling.  The  reverse  is  true  for  thick 
ice.  However,  in  both  these  cases  the  three-dimensional  stress  will  increase  as  the  ice  becomes 
thicker.  This  is  in  contrast  to  the  rubble  case  where  the  three-dimensional  stress  will  be  a  con¬ 
stant  at  *  1.6x  10*  N  m'J  (or  *  2.3  lb  in.'J)- 


NUMERICAL  SIMULATION  RESULTS 


The  previous  section  has  analytically  examined  the  character  of  the  ridge  redistribution 
process.  However,  to  examine  the  performance  of  this  assumed  process  when  it  is  coupled  to 
ice  dynamic  and  thermodynamic  equations,  it  is  necessary  to  carry  out  numerical  Simula- 


ridge  made  from  a  single  thickness  of  ice  would  have  vertical  sides.  The  Hibler  redistributor 
avoids  this  feature  by  uniformly  distributing  ice  from  twice  the  block  thickness  up  to  some 
maximum  thickness.  Also,  perhaps  more  importantly,  the  maximum  thickness  is  taken  to 
scale  with  the  square  root  of  the  thickness  of  ice  being  ridged.  This  idea  reflects  the  physical 
notion  that  for  equal  amounts  of  deformation  per  ridge,  doubling  the  ice  thickness  will  not 
double  the  maximum  ridge  height.  This  type  of  scaling  is  supported  by  observations  (Tucker 
and  Govoni  1981).  The  modified  Hibler  case  is  a  uniform  redistributor  without  this  scaling. 
The  rubble  redistribution  is  an  idealized  case  in  which  the  same  thickness  of  ice  is  always 
created. 

Comparison  to  ridge  morphological  data 

The  scaling  characteristics  of  different  redistributors  can  be  tested  by  field  examination  of 
ridge  height  and  block  size  characteristics.  A  particularly  useful  data  set  for  this  purpose  has 
been  obtained  and  analyzed  by  Tucker  and  Govoni  (1981).  This  data  set  was  taken  off  the 
North  Slope  of  Alaska.  Figure  3  shows  the  salient  characteristics  of  the  data  set  that  are  rele¬ 
vant  to  the  redistribution  theory.  The  error  bars  are  the  standard  error  in  the  mean  ridge 
height  estimate.  As  can  be  seen  from  this  figure,  the  square  root  scaling  fits  the  data  best. 

To  approximately  convert  the  ridge  height  scaling  to  an  ice  thickness  scaling,  a  4: 1  ridge 
keel  to  sail  ratio  is  assumed  (Kovacs  and  Mellor  1974).  With  this  scaling  the  square  root  case 
is  equivalent  to  an  H*  of  100  m  in  the  Hibler  redistributor,  whereas  the  linear  scaling  repre¬ 
sents  a  factor  of  K  =  15  in  the  Thorndike  et  al.  or  modified  Hibler  redistributors.  The  rubble 
case  represents  a  fixed  thickness  H'  =  20  m. 
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Figure  3.  Comparison  of  different  redistribution  scaling  laws  with  observed 
data.  The  solid  curve  represents  the  Hibler  scaling,  the  long-dashed  curve  the 
Thorndike  et  al.  or  modified  Hibler  scaling,  and  the  short -dashed  curve  the  rubble 
scaling. 

Ice  strength  for  different  redistribntors 

Of  particular  relevance  to  ice  mechanics  are  the  ice  strengths  generated  by  different  redis¬ 
tributors.  To  get  some  feeling  for  these  strengths,  consider  the  special  case  of  only  one  thick¬ 
ness  of  ice,  say  A«,  being  ridged.  Formally  representing  this  condition  by  F\h)  -  S(h  -  h,), 
one  can  substitute  eq  16  into  eq  18  and  obtain  analytical  results  for  the  ice  strength.  After 
some  algebra,  different  redistributors  yield  the  following  strength  equations: 


where  g(h)dh  is  the  fraction  of  area  covered  by  ice  with  thickness  between  h  and  h  +  dh.  If  we 
consider  a  pure  convergent  deformation  with  convergence  rate  V  •  £  ,  then  -  V  .  v  Wr(h,g) 
is  the  change  in  the  ice  thickness  distribution  g(h)  per  unit  time  because  of  deformation. 

In  eq  16,  P(h)  is  a  probability  function  specifying  which  categories  are  being  destroyed  by 
ridging.  The  quantity  y(h  ',h)dh  can  be  thought  of  as  the  area  of  ice  put  into  thickness  inter¬ 
val  (h,  h  +  dh)  when  a  unit  area  of  ice  thickness  h '  is  used  up.  Basically,  the  integral  over  y 
specifies  where  the  ice  is  transferred  during  ridging  by  changing  the  areal  thickness  distribu¬ 
tion  function  g(h).  It  is  this  redistribution  function  y  that  is  of  particular  relevance  to  this 
paper.  It  can  be  shown  (Appendix  A)  that  to  conserve  mass,  y(h \h)  must  obey  the  equation 
00 

J*  y(h',h)hdh  =  h‘.  (17) 

0 

Strengths  can  be  estimated  from  this  kind  of  theory,  if  it  is  insisted  that  the  rate  of  defor¬ 
mation  work  equals  the  work  done  through  ridge  building  (Rothrock  1975).  This  constraint 
(see  Appendix  A)  leads  to  the  equation  for  the  two-dimensional  ice  strength  P *: 

00 

P*  =  C  f  W,  hl  dh  (18) 

0 

where  the  constant  C 
C  =  c 

e» 

is  obtained  by  assuming  the  change  in  gravitational  potential  energy  during  ridging  is  equal 
to  the  frictional  losses  in  ridging.  In  this  expression  q,  and  g.  are  ice  and  water  densities  and 
g  is  the  acceleration  of  gravity. 

Some  specific  redistribntors 

Several  different  forms  for  y(huh,)  have  been  suggested.  Others,  although  not  previously 
suggested,  are  useful  for  pedagogical  reasons.  A  variety  of  redistributors  satisfying  eq  17  are 
listed  below. 

y(h„hi)  =  yih.-Kh,)  (UK) 
with  K  =  constant  (Thorndike  et  al.), 

r(A.,A.)  =  (l/2(//*-A,)J  for  2h{  sh,s2  yfH*  -JK, 
with  H*  =  constant  (Hibler), 

y(h„hi)  =  [2//i,(Ar2 -4)]  for  2ht  <,  h,  <,  Kh , 
with  K  =  constant  (modified  Hibler), 
y(huh>)  =  m-H')(h,/H') 

(rubble  redistribution). 

The  Thorndike  et  al.  (1975)  redistributor  specifies  that  ice  is  transferred  into  a  fixed  multi¬ 
ple  of  its  thickness.  While  this  is  computationally  simple  it  has  the  unrealistic  feature  that  a 


or  ridging  ice  if 


oo 

f  gdh  >  1. 

o 

(Note  that  while  g  is  normalized  at  the  beginning  of  the  time  step,  it  will  not  be  after  eq  1 1 .) 
As  part  of  this  redistribution  step,  any  negative  g  values  (because  of  second  order  differenc¬ 
ing  in  the  horizontal  advection  equations)  are  set  equal  to  zero.  The  net  heat  imbalance  re¬ 
sulting  from  this  removal  is  used  to  increase  the  mixed  layer  temperature  to  ensure  that  the 
heat  budget  is  conserved.  The  final  values  of  mixed  layer  temperature  and  g  are  obtained  in 
eq  14  and  15  by  laterally  melting  ice  until  either  no  ice  is  left  or  the  mixed  layer  temperature  is 
at  freezing.  Since  any  ice  area  lost  through  melting  will  be  compensated  for  by  increasing  the 
open  water  fraction,  these  steps  ensure  that  g  is  normalized  at  the  end  of  the  time  step.  Note 
that  in  this  sequence  deformation  effects  are  included  by  the  $,(u/i)  terms  and  by  the  open 
water  creation.  These  in  turn  can  necessitate  ridging.  Handling  the  deformation  in  this  some¬ 
what  indirect  manner  allows  generalization  to  an  Eulerian  grid. 


ANALYTIC  EXAMINATION  OF  THE 
RIDGE  REDISTRIBUTION  PROCESS 

On  the  geophysical  scale  stresses  in  pack  ice  are  largely  determined  by  the  ridge  building 
process.  To  model  this  process  on  a  large  scale  it  has  been  suggested  that  thin  ice  be  redis¬ 
tributed  into  thicker  ice  categories  (see  Appendix  A).  The  precise  manner  in  which  this  redis¬ 
tribution  should  be  carried  out  is  not  clear.  As  an  initial  guess,  Thorndike  et  al.  (1975)  sug¬ 
gested  a  redistribution  that  transfers  ice  into  categories  that  are  a  fixed  multiple  of  the  initial 
thickness.  Ridge  observations  and  theoretical  considerations  suggest  that  such  a  redistribu¬ 
tion  is  unrealistic.  In  order  to  provide  a  more  realistic  parameterization  of  the  ridging  pro¬ 
cess,  a  scaling  law  for  ridge  building  is  proposed  here  (see  Appendix  A).  The  purpose  of  this 
section  is  to  examine  in  more  detail  the  strength  characteristics  of  different  redistribution 
processes  (including  the  particular  scaling  law  proposed  here),  and  to  compare  these  process¬ 
es  to  ridge  morphological  data. 

Theoretical  framework 

In  variable  thickness  sea  ice  models,  the  ridging  process  is  parameterized  by  probability 
functions  specifying  how  different  ice  categories  are  modified  by  deformation.  Following  the 
notation  used  in  Appendix  A,  the  function  lVr(h,g)  describes  this  ridging  process: 

00  00 

W*h,g)  =  l-P(h)  g(h)+  J  y  (h  \h )  P{h ')  g(h ')  dh')/J  [ P(h)  g(h )  - 

0  0 

oo 

/  y(h',h)Pih')g(h')dh']dh 
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(16) 


the  redistributor  that  transfers  thin  ice  to  thicker  categories,  and  Tis  the  boundary  layer  tem¬ 
perature.  Also 


/<*)  =  /  Ah)  g(h)  dh 

0 

oo 

/*(g)  =/  Mh)  g(h)  dh 
0 

oo 

FL(g,  T)  =  J*  FL(g,  h,  T)  dh. 

0 

For  illustration,  only  the  x  component  of  the  ice  velocity  field  (denoted  by  u)  has  been  used 
and  the  freezing  point  of  seawater  has  been  set  equal  to  zero.  Also,  the  Du  term  schematical¬ 
ly  represents  a  horizontal  diffusion  term.  In  the  complete  code  these  diffusion  terms  are 
identical  to  those  used  in  the  two-level  model  (Hibler  1979)  and  were  included  for  suppres¬ 
sion  of  nonlinear  instabilities.  In  time,  u  is  defined  at  (*  “  and  the  other  variables  in  eq  8  and 
9  at  t‘.  At  the  beginning  of  the  time  step  it  is  assumed  that  the  vertical  growth  rates  are  ob¬ 
tained  from  eq  5.  To  advance  g  and  T requires  the  following  steps  (spatial  and  thickness  dif¬ 
ferences  are  denoted  by  5,  and  6„,  time  locations  by  superscripts,  and  intermediate  values 
within  the  time  step  by  subscripts): 


=  *‘-A/5W) 

00) 

-  *'  ~  A  t («,[«' + *(«'  +  g‘, )0.5]  +  +  Dg‘  -  Mu' " ")  1 

(ID 

-  r+Af(Ag‘) -/»(*')! 

(12) 

=  g!+A/*,  (*J+,,«0 

(13) 

=  g)  + At  FL(g'*l,T!) 

(14) 

=  Ti  +  AtFL(g{",  TO. 

(15) 

In  this  sequence  eq  10  provides  a  provisional  value  of  g  to  approximately  center  the  spatial 
advection  term.  This  approach  supplies  a  modified  Euler  step  (Kurihara  1965),  which  is  sec¬ 
ond  order  accurate  in  time  for  the  horizontal  advection.  After  this  takes  place  the  remaining 
steps  are  basically  sequential  splitting  steps.  Equation  1 1  determines  the  change  of  g  due  to 
spatial  advection,  thickness  advection  and  changes  in  open  water.  All  advection  terms  are 
done  conservatively  so  that  all  g  values  are  conserved  over  the  global  grid.  In  addition,  as  dis¬ 
cussed  in  Appendix  C,  the  thickness  advection  term  is  done  in  a  manner  that  conserves  the 
surface  heat  budget  components  and  is  stable  for  forward  time  steps.  In  eq  12  any  heat  ab¬ 
sorbed  by  open  water  that  is  not  used  in  the  vertical  growth  rates  is  used  to  increase  the  mixed 
layer  temperature.  In  eq  13  the  mechanical  redistribution  process  is  carried  out  (see  Appen¬ 
dix  A).  This  consists  of  normalizing  g  by  either  creating  more  open  water  if 
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Figure  2.  Schematic  arrangement  of  thickness  partition  used  in  numerical  calculations. 

growth  by  Martin  (1979)  shows  that  growth  of  thin  ice  can  be  substantially  enhanced  by 
frazil  ice  formation.  Such  ice  is  derived  from  small  ice  crystals  that  form  in  the  water  col¬ 
umn.  Recent  observations  by  Ackley  et  al.  (1980)  of  the  Antarctic  sea  ice  field  indicate  that 
the  frazil  ice  component  of  Weddell  Sea  pack  ice  is  substantially  larger  than  that  found  in  the 
Arctic.  As  a  consequence,  inclusion  of  enhanced  thin  ice  growth  may  be  important  for  Ant¬ 
arctic  pack  ice  simulations. 

Numerical  scheme 

Equations  4-7  are  numerically  solved  as  an  initial  value  problem  using  finite  difference 
techniques.  While  the  strengths  are  generated  differently,  the  numerical  coupling  of  these 
equations  with  the  dynamics  proceeds  in  the  same  way  as  employed  previously  for  a  coupled 
two-level  model  (Hibler  1979).  Since  the  numerical  scheme  for  the  two-level  model  is  dis¬ 
cussed  by  Hibler  (1979),  it  is  enough  here  to  present  the  time  marching  scheme  and  finite  dif¬ 
ferences  used  to  solve  eq  4-7.  The  time  marching  scheme  is  briefly  sketched  below.  Details 
on  the  thickness  partition  and  finite  difference  code  are  given  in  Appendix  C.  In  general  the 
numerical  complexity  of  the  equations  is  substantially  increased  by  the  inclusion  of  the  ad- 
vection  terms.  These  advection  terms  arise  because  fixed  grids  in  both  space  and  thickness 
are  used.  The  spatial  grid  is  a  regular  mesh  staggered  configuration  (see  Fig.  5,  Hibler  1979) 
with  thickness  characteristics  defined  at  the  center  and  velocity  points  defined  at  the  corners. 
The  thickness  grid  is  composed  of  an  arbitrary  number  of  irregularly  spaced  thickness  levels. 
For  the  simulations  performed  here,  ten  thickness  levels  are  employed  as  shown  in  Figure  2. 
In  this  grid  larger  categories  are  used  for  thicker  ice  since  growth  rates  change  less  rapidly 
there.  Open  water  is  considered  simply  by  having  one  thickness  category  centered  at  zero 
thickness.  These  thickness  levels  remain  fixed.  However,  the  relative  areal  extent  of  the  ice  in 
each  category  evolves  in  response  to  vertical  ice  growth,  spatial  advection,  deformation  and 
lateral  melting. 

To  solve  eq  4-7  a  series  of  splitting  operations  is  used.  In  a  fixed  Eulerian  grid,  the  se¬ 
quence  in  which  these  operations  are  performed  is  critical  since  it  is  important  that  certain 
operations  be  done  last  to  ensure  conservation  of  several  quantities.  The  time  marching 
scheme  for  the  coupled  thickness  equations  can  be  illustrated  using  the  following  simplified 
equations: 

O/a  0=  -(d/dx)(ug)-(d/dh)(fg)  +  Fl(g,T)  +  t,(g,u)  +  t,(u)-Du  (8) 

0  T/d  t)  =  FL(g,T)  +  [f  {g)~  fb(g)]  (9) 

where  is  the  portion  of  the  redistributor  that  only  creates  open  water,  ,  is  the  portion  of 
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tion  to  the  amounts  of  heat  absorbed  by  leads  and  not  used  in  increasing  the  rate  of  vertical 
melt.  Equation  7,  on  the  other  hand,  states  that  there  will  be  enough  instantaneous  lateral 
melt  to  either  lower  the  mixed  layer  temperature  to  freezing  or  to  remove  all  the  ice.  Also,  eq 
7  assumes  that  lateral  melt  reduces  all  ice  thicknesses  by  the  same  percentage.  The  physical 
argument  here  is  that  ice  will  lose  mass  by  lateral  melting  in  proportion  to  its  thickness.  This, 
in  turn,  is  based  on  the  physical  notion  that  thicker  ice  will  have  a  larger  vertical  interface 
with  the  ocean  than  thinner  ice.  Coupling  the  equations  in  this  way  removes  the  necessity  of 
considering  any  special  cases.  For  example,  cooling  off  a  warm,  ice-free  mixed  layer  is  natu¬ 
rally  treated  in  eq  5-7  by  simultaneously  freezing  ice  and  removing  the  frozen  ice  by  lateral 
melting. 

While  better  parameterizations  are  needed,  numerical  experiments  with  this  mixed  layer 
growth  model,  coupled  to  the  variable  ice  thickness  model  discussed  earlier,  can  help  deter¬ 
mine  the  relative  role  of  dynamics  and  thermodynamics  in  sea  ice  growth,  drift  and  decay. 
The  idea  in  the  (0)  aO  case  is  that  under  “winter  conditions”  the  growth  rates  of  the  sea  ice 
cover  are  basically  vertical  and  are  dominated  by  heat  budget  considerations.  This  assump¬ 
tion  is  reasonable  in  the  Arctic  Basin,  where  leads  do  not  long  remain  ice-free,  and  is  consist¬ 
ent  with  sensitivity  studies  by  Maykut  (1978).  The  heat  budget  code  described  here  approxi¬ 
mates  these  features.  The  main  shortcoming  is  that  the  effect  of  snow  cover  on  the  conduc¬ 
tivity  of  the  snow-ice  system  is  not  included.  While  this  feature  is  not  critical  to  equilibrium 
thicknesses  (see,  e.g.,  Maykut  and  Untersteiner  1971),  snow  does  substantially  affect  the 
growth  rates  of  thin  ice  (Maykut  1978).  Improvements  of  this  parameterization  are  needed. 
However,  it  should  be  noted  that  in  a  coupled  model  inadequate  estimates  of  thin  ice  growth 
rates  may  well  be  offset  by  larger  thin  ice  growth  estimates  shortening  the  thin  ice  lifetime. 
Consequently,  time-averaged  energy  exchange  may  not  be  critically  affected. 

Under  melting  conditions  the  open  water  heat  absorption  terms  reflect  the  fact  that 
boundary  layer  warming  and  lateral  melting  can  substantially  affect  the  decay  of  sea  ice.  Ef¬ 
fects  of  this  kind  have  been  discussed  in  reviews  by  Rothrock  (1979),  Wadhams  (1980)  and 
Hibler  (1980b).  Lateral  melting  has  been  particularly  emphasized  by  Zubov  (1945)  and  Lan- 
gleben  (1972)  in  studies  of  shore-fast  ice  decay.  These  authors  have  explained  observation  of 
shore-fast  ice  decay  by  using  absorbed  radiation  solely  in  decreasing  the  horizontal  dimen¬ 
sions  of  floes.  Wadhams  (1980)  has  suggested  that  near  the  ice  edge,  lateral  melting  may  de¬ 
pend  upon  the  geometric  properties  of  the  ice.  Specifically,  as  floes  become  smaller  and  more 
numerous,  the  effective  “length”  of  flow  edges  increases.  In  addition  to  lateral  melting,  the 
movement  of  the  ice  and  the  wind  mixing  will  cause  the  ice  to  move  over  regions  that  were 
formerly  leads.  Summer  observations  of  the  mixed  layer  and  ice  drift  by  McPhee  (1980b)  in¬ 
dicate  that  substantial  vertical  mixing  can  be  induced  by  the  wind  and  ice  motion.  Also,  near 
the  ice  edge  turbulence  and  wave  effects  can  cause  horizontal  mixing  and  enhance  vertical 
melt  (see,  e.g.,  Wadhams  et  al.  1979).  The  parameterization  suggested  here  approximates 
these  processes  primarily  by  increasing  vertical  melting  due  to  heat  absorption.  In  addition, 
for  a  sufficiently  loose  pack,  lateral  melting  terms  are  added. 

In  both  the  decay  and  growth  cases  a  possible  modification  would  be  to  allow  different 
growth  and  decay  rates  for  thick  ice.  Very  thick  ice  in  the  form  of  ridges  may  well  behave  dif¬ 
ferently  than  level  ice.  Ablation  observations  by  Koerner  (1973),  for  example,  indicate  that 
the  upper  surfaces  of  first-year  ridges  ablate  much  more  rapidly  than  level  ice.  Such  effects 
may  be  even  more  pronounced  on  the  bottoms  of  ridges,  where  the  deep  keel  of  a  ridge  al¬ 
lows  ablation  at  the  sides  as  well  as  on  the  bottom  of  the  ice.  In  light  of  these  considerations 
some  type  of  parameterization  of  the  growth  rates  as  a  function  of  the  geometric  properties 
of  pack  ice  is  needed. 

Finally,  in  the  growth  case  it  is  possible  that  thin  ice  growth  rates  may  be  substantially 
larger  than  those  based  on  heat  budget  considerations.  Specifically,  analysis  of  young  ice 
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Under  conditions  where  open  water  is  losing  heat  to  the  atmosphere,  the  heat  budget 
growth  rates  [denoted  by  Mh)\  are  taken  to  be  the  vertical  growth  rates  in  eq  4.  However, 
when  open  water  is  absorbing  heat,  the  heat  is  allowed  to  mix  underneath  the  ice  flow  and 
decrease  the  vertical  growth  rate  of  all  categories  down  to  some  minimum  value.  Any  re¬ 
maining  heat  is  allowed  to  either  cause  lateral  melting  or  raise  the  temperature  of  the  mixed 
layer.  In  this  specific  parameterization  the  mixed  layer  temperature  is  always  kept  at  freezing 
in  the  presence  of  an  ice  cover.  Consequently,  all  available  heat  absorbed  by  leads  and  not 
used  in  vertical  melting  is  used  in  lateral  melting  until  the  ice  disappears.  Also,  under  growth 
conditions  no  ice  is  allowed  to  form  until  the  mixed  layer  reaches  the  freezing  temperatures 
of  seawater. 

The  vertical  rates  in  this  parameterization  may  be  formally  described  by 


m  » 


Max[0,  /»(<))] 

Max|/,, [/.(A)  +  Min(/6(0),0)J| 


if  ft  =  0 


if  h  >  0 


where  g,  is  the  fraction  of  area  covered  by  ice,  which  is  defined  as 


g,  =  Lim  f  g(h)dh. 

t-0  J 


Specifying  a  minimum  decay  rate,  /,,  is  based  on  the  physical  notion  that  there  is  a  limit  to 
how  rapidly  sea  ice  may  melt  in  water  near  the  freezing  point.  In  the  simulation  performed 
here,/,  is  taken  to  be  -15  cm/day.  This  number  is  close  to  the  maximum  wave-induced  melt¬ 
ing  rates  calculated  by  Wadhams  et  al.  (1979). 

The  lateral  melting  and  warming  of  the  muted  layer  are  described  by  the  coupled  equations 

=  m.2  +  [QACjdmi')\\]  FL{g,h,Tml.)hdh  +  /  W)~Mh)\g(h)dl\  (6) 


FL(h,g,Tm„)  = 


~C{Tml,)g(h) 


CiTm„)gM.h) 


for  h  >0 


for  A  =  0 


where 


C(Tml.)  =  Mini— Cjd„u,  I 

I 

with  Ti  the  mean  ice  thickness 


h  =  J*  g(h)hdh. 


In  these  equations  Tm„  is  the  mixed  'ayer  temperature  in  Kelvins,  dmU  is  the  mixed  layer 
depth,  Q,  is  the  volumetric  heat  of  fusion  of  ice  (set  at  302  MJ  m*’),  and  C.  is  the  volumetric 
heat  capacity  of  water  (set  at  4.19  MJ  m'1)-  Also,  in  eq  6  and  7  the  freezing  temperature  of 
seawater  is  taken  to  be  271.2  K.  Equation  6  specifies  that  the  mixed  layer  temperature  will 
decrease  in  proportion  to  the  mass  of  ice  removed  by  lateral  melting  and  increase  in  propor- 
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causes  the  more  general  thickness  distribution  equation  to  become 


is 

dt 


+  V.(tig)  + 


d(fg) 

dh 


Fl  + 


(4) 


Considerable  complexity  in  the  equation  resides  in  the  mechanical  redistribution  function 
\p.  Because  of  this  complexity  a  more  complete  discussion  of  ^  is  given  in  Appendix  A. 
There,  a  specific  redistribution  function  is  proposed  and  the  relation  of  the  ice  strength  to  the 
ice  thickness  distribution  is  briefly  discussed.  In  addition,  because  of  the  importance  of  ridg¬ 
ing,  the  ramifications  of  the  ridge  distribution  process  proposed  in  Appendix  A  are  analyti¬ 
cally  examined  and  compared  to  observed  data  in  a  later  section.  However,  for  the  present 
discussion  it  is  enough  to  note  that  \p  creates  open  water  under  pure  divergence  and  transfers 
thin  ice  to  thick  ice  categories  under  pure  convergence.  For  an  arbitrary  deformation  state, 
the  open  water  creation  and  ridging  are  considered  to  be  simultaneous.  In  addition  the 
strength  of  the  system  is  related  to  the  amount  of  thin  ice,  so  that  thicker,  more  compact  ice 
will  exhibit  greater  strength  and,  hence,  greater  resistance  to  compression. 

Heat  budget  and  oceanic  boundary  layer 

Given  vertical  and  lateral  growth  rates,  one  could  use  the  thickness  distribution  (eq  4)  to 
determine  the  evolution  of  g.  However,  to  determine  these  growth  rates  is  a  complex  task 
and  it  minimally  requires  consideration  of  the  oceanic  boundary  layer  and  heat  budget  at  the 
air/sea  ice  interface.  For  this  purpose  a  simple  parameterization  consistent  with  available 
knowledge  is  suggested  here.  Some  of  the  complexity  that  might  be  considered  in  a  more  de¬ 
tailed  parameterization  is  discussed  later  in  this  section. 

Basically,  at  each  ice  thickness  category  vertical  growth  rates  are  estimated  1)  from  heat 
budget  considerations  at  the  top  and  bottom  surface  of  the  ice  and  2)  by  adding,  via  lateral 
mixing,  some  of  the  heat  absorbed  by  leads.  In  the  heat  budget  calculations  a  simple  time- 
independent  thermodynamic  sea  ice  model  is  used.  Following  Semtner  (1976),  this  model  ap¬ 
proximates  the  heat  transfer  through  the  ice  by  assuming  a  linear  temperature  profile  to¬ 
gether  with  a  constant  ice  conductivity.  However,  in  contrast  to  Semtner  (1976),  the  conduc¬ 
tivity  effects  of  a  snow  cover  are  not  considered.  Instead,  following  Bryan  et  al.  (1975)  and 
Manabe  et  al.  (1979)  the  model  approximates  the  effects  of  snow  cover  by  allowing  the  ice 
surface  albedo  to  be  that  of  snow  when  the  calculated  surface  temperature  is  below  freezing 
and  that  of  snow-free  ice  wHn  the  surface  temperature  is  at  the  melting  point.  With  these  as¬ 
sumptions  the  upward  heat  flux  (/*)  through  ice  of  thickness  H  is 

/„  =  (K/H)  (71  -  r„) 

where  K  is  the  ice  conductivity,  Tw  the  water  temperature  and  T.  the  surface  temperature  of 
the  ice.  This  simple  thermodynamic  model  is  used  in  conjunction  with  a  surface  heat  budget 
computation  similar  to  that  of  Parkinson  and  Washington  (1979)  and  Manabe  et  al.  (1979). 
In  this  computation  (see  Appendix  A)  the  surface  temperature  of  the  ice  that  balances  the 
surface  heat  budget  is  obtained  by  iteration.  This  temperature  then  dictates  the  conduction 
of  heat  through  the  ice  and,  hence,  the  ice  growth  rate.  If  the  iteration  yields  an  above- 
freezing  value,  the  surface  temperature  is  set  at  the  freezing  point.  Surface  and  bottom  abla¬ 
tion  rates  are  then  determined  by  imbalances  in  the  surface  heat  budget  and  by  the  conduc¬ 
tion  of  heat  into  the  mixed  layer.  As  in  Maykut  and  Untersteiner  (1971),  the  effect  of  heat 
transfer  from  deep,  warmer  ocean  water  is  approximated  by  assuming  a  constant  oceanic 
heat  flux  into  the  mixed  layer.  This  thermodynamic  model  and  surface  heat  budget  computa¬ 
tion  is  described  in  more  detail  in  Appendix  B. 


ice  rheology  a  viscous  plastic  constitutive  law  is  used.  Rigid  plastic  behavior  is  approximated 
in  this  law  by  allowing  the  ice  to  flow  plastically  for  normal  strain  rates  and  to  creep  in  a 
linear  viscous  manner  for  small  strain  rates.  Documentation  of  this  dynamical  code  is  avail¬ 
able  (Hibler  1980a). 

In  the  work  presented  here,  the  dynamical  formulation  used  by  Hibler  (1979)  is  employed. 
However,  a  more  general  treatment  of  the  ice  thickness  distribution,  the  ice  strength  and  the 
thermodynamic  code  is  presented. 

Ice  thickness  equations 

For  the  purposes  of  this  paper  the  ice  is  considered  to  be  a  two-dimensional  continuum 
with  a  velocity  field  u.  This  velocity  field  is  obtained  by  solving  a  momentum  balance  (see 
Hibler  1979)  that  includes  the  effects  of  internal  ice  stress.  The  stress  state  in  the  ice  is  de¬ 
scribed  by  a  plastic  constitutive  law 

<*ij  =  (1) 

where  is  the  two-dimensional  stress  tensor,  e0  the  strain  rate  tensor  and  P*  the  ice 
strength.  As  discussed  in  more  detail  in  Appendix  A,  P*  is  a  function  of  the  amount  of  thin 
ice  present  in  the  ice  cover. 

As  in  Thorndike  et  al.  (1975),  the  ice  thickness  characteristics  are  described  by  an  areal  ice 
thickness  distribution  g(h,£,t)  where  g(h,)c,t)dh  is  the  fraction  of  area  (in  a  region  centered 
at  position  £  at  time  t)  covered  by  ice  with  thickness  between  h  and  h  +  dh.  This  distribution 
evolves  in  response  to  deformation,  advection,  growth  and  decay.  Neglecting  lateral  melting 
effects,  Thorndike  et  al.  (1975)  derived  the  following  governing  equation  for  the  thickness 
distribution: 

37  +  *•<«*>  +TT~*  ® 

where  /is  the  vertical  growth  (or  decay)  rate  of  ice  of  thickness  h  and  ^  is  a  redistribution 
function  (depending  on  h  and  g)  that  describes  the  creation  of  open  water  and  the  transfer  of 
ice  from  one  thickness  to  another  by  rafting  and  ridging.  In  general,  /is  a  function  of  g(h), 
and  $  is  a  function  of  both  g(h)  and  the  rate  of  deformation.  Except  for  the  last  two  terms, 
eq  2  is  a  normal  continuity  equation  for  g.  The  last  term  on  the  left-hand  side  can  also  be 
considered  a  continuity  requirement  in  thickness  space  since  it  represents  a  transfer  of  ice 
from  one  thickness  category  to  another  by  the  growth  rates.  An  important  feature  of  the 
Thorndike  et  al.  (1975)  theory  is  that  it  presents  an  “Eulerian”  description  in  thickness 
space.  In  particular,  growth  takes  place  by  rearrangement  of  the  relative  areal  magnitudes  of 
different  thickness  categories. 

Equation  2  can  be  formally  generalized  to  include  lateral  growth  by  simply  adding  a  source 
and  sink  term  FL(^,h)  such  that 


00 

f  FLdh  =  0.  (3) 

o 

Equation  3  follows  from  ti.e  fact  that,  by  definition,  lateral  melting  (or  freezing)  of  ice  will 
be  compensated  for  by  a  change  in  the  extent  of  open  water.  In  addition  FL  >  0  for  h  =  0 
and  Fl  <  0  for  h  >  0.  These  conditions  reflect  lateral  melting  decreasing  the  areal  extent  of 
ice  while  increasing  the  open  water  extent.  The  addition  of  this  lateral  growth  term  to  eq  2 
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Figure  7  (cont’d).  Average  April  and  August  thickness  con¬ 
tours  (m)  for  the  fifth  year  of  the  standard  and  ‘‘thermo¬ 
dynamics-only"  simulations.  The  dashed  line  represents  the 
average  40%  ice  concentration  contour  at  the  end  of  August  for  the 
time  period  1966-1976  as  estimated  by  the  British  Meteorological 
Office.  Bracknell,  England. 


fects  are  very  similar  to  those  obtained  by  Washington  et  al.  (1976),  with  ice  exceeding  3  m  in 
the  central  basin  in  winter,  and  all  the  ice  at  certain  near-shore  locations  melting  in  summer. 
The  standard  ca$e,  on  the  other  hand,  exhibits  a  pronounced  ice  buildup  along  the  Canadian 
Archipelago.  The  reason  for  the  spatial  variation  in  the  standard  case  is  apparent  from  ex¬ 
amination  of  the  annual  ice  velocity  field  in  Figure  8.  This  field  exhibits  a  clockwise  Beaufort 
Gyre  and  transpolar  drift  in  agreement  with  observations  of  the  mean  annual  drift  (see,  e.g., 
Gordienko  1958).  While  instantaneous  velocities  may  differ  substantially,  on  the  average  the 
drift  tends  to  thin  the  ice  off  the  Alaskan  and  Siberian  coasts  while  increasing  the  ice  thick¬ 
ness  off  the  Canadian  Archipelago.  The  shape  and  magnitude  of  the  standard  case  ice  thick¬ 
ness  contours  agree  well  with  observed  estimates.  Figure  9,  for  example,  shows  observed 
results  based  on  submarine  data,  portions  of  which  are  reported  by  LeShack  et  al.  (1971). 
Similar  spatial  thickness  variations  have  been  obtained  from  airborne  and  subsurface  studies 
of  the  ridge  statistics  in  the  Arctic  Basin  (see,  e.g.,  Hibler  et  al.  [1974]  and  Wadhams  [1981]). 

The  summer  ice  edge  characteristics  of  the  standard  case  are,  however,  less  realistic.  Spe¬ 
cifically,  while  the  shape  of  the  ice  edge  is  reasonable,  its  location  is  consistently  too  far  from 
shore.  It  also  appears  that  the  standard  case  is  yielding  rather  large  amounts  of  open  water  in 
the  central  basin.  This  feature  is  demonstrated  in  Figure  10,  which  shows  mid-month  com¬ 
pactness  transects  from  July  through  October.  These  transects  were  taken  along  a  line  con¬ 
necting  grid  cells  4  and  8  (see  Fig.  5).  Using  a  large  amount  of  data  from  aerial  reconnais¬ 
sance  flights,  Wittmann  and  Schule  (1966)  have  estimated  that  from  August  to  October  there 
is  typically  about  12%  open  water  in  the  Canadian  Basin  (between  the  Pole  and  Alaska),  and 
7%  in  the  Eurasian  Basin  (between  the  Pole  and  Franz  Josef  Land).  The  August  through  Oc¬ 
tober  averages  of  the  simulated  values  agree  with  the  geographical  trend  of  these  observa¬ 
tions,  but  are  somewhat  larger.  However,  the  substantial  month-to-month  variations  in  the 
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Figure  8.  Average  annual  ice  velocity  field  for  the  fifth  year  of  the 
standard  simulation.  A  velocity  vector  of  one  grid  space  long  represents 
0.02  msr\ 
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Figure  9.  Approximate  contours  ( m )  of  observed  ice 
thickness  values  obtained  from  submarine  sonar  data 
(personal  communication  with  L.A.  LeShack,  LeShack 
Associates,  Silver  Spring,  Maryland,  1979).  The  solid 
contours  were  obtained  by  LeShack  from  a  composite  analysis 
of  both  summer  (I960,  1962)  and  winter  (I960)  data,  whereas 
the  dashed  contours  are  from  April  1977  data. 
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Figure  10.  Mid-month  ice  compactness  transects 
for  the  fifth  year  of  the  standard  simulation. 
These  transects  were  taken  along  a  line  connecting  grid 
cells  4  and  8. 


simulated  compactness  values  (see  Fig.  10)  indicate  that  real  time  observations  are  needed  for 
detailed  verification.  Similar  considerations  also  apply  to  ice  edge  variations,  since  year-to- 
year  values  of  the  ice  edge  can  differ  significantly  from  climatological  means  (see,  e.g., 
Walsh  and  Johnson  1979b). 

Of  particular  interest  to  this  author  are  the  ridge  buildup  results.  While  the  model  does  not 
explicitly  keep  track  of  ridging,  it  is  possible  to  determine  the  volume  of  deformed  ice  created 
over  an  annual  cycle  at  a  fixed  location.  Figure  1 1  shows  contours  of  this  quantity.  The  gen¬ 
eral  shape  of  the  roughness  contours  agrees  well  with  surface  roughness  observations  in  the 
western  Arctic  Basin  (Weeks  et  al.  1971,  Hibler  et  al.  1974).  The  two  major  features  are  a 
heavy  buildup  of  ridging  off  the  Canadian  Archipelago  and  less  ridging  in  the  Beaufort  Sea 
than  near  the  Pole.  Note,  however,  that  there  is  a  zone  of  heavy  ridging  just  off  the  North 
Slope,  in  agreement  with  observations  by  Tucker  et  al.  (1979).  The  other  major  feature  is  a 
tongue  of  high  ridging  further  offshore  near  the  tip  of  Greenland.  Roughness  data  obtained 
from  submarine  sonar  profiles  by  LeShack*  also  show  such  a  tongue. 

Ice  edge  evolution  and  sensitivity 

The  ice  edge  characteristics  take  several  years  to  fully  evolve.  This  evolution  time  is  illus¬ 
trated  in  Figure  12,  which  shows  the  mid-September  ice  thicknesses  along  a  transect  perpen¬ 
dicular  to  the  North  Slope.  As  is  apparent,  the  ice  edge  takes  about  3  years  to  completely 
develop,  with  the  difference  between  the  first  two  summers  being  most  pronounced.  Because 
of  this  long  evolution  time,  it  is  possible  that  some  of  the  unrealistic  results  may  be  attribut¬ 
able  to  the  choice  of  the  wind  field  for  a  particular  year.  By  considering  interannual  variabil¬ 
ity  in  a  simulation,  earlier  years  with,  say,  less  persistent  offshore  winds  in  summer  could  re¬ 
duce  the  formation  of  an  excessive  ice  edge  in  a  particular  year. 


Ontario*  from  Coo«*  (km) 


Figure  12.  Evolution  of  the  mid-September 
( day  257)  ice  thickness  characteristics  in  the 
standard  simulation.  Thickness  transects  along  a 
line  connecting  grid  cells  4  and  8  are  plotted. 


Personal  communication  with  L.A.  LeShack.  LeShack  Associates,  Silver  Springs,  Maryland,  1979. 
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Figure  13.  Average  August  thickness  contours  (m)  for  the 
high-strength  and  high-growth  simulations. 


However,  while  the  interannua!  variability  may  be  a  factor,  it  seems  likely  that  the  unreal¬ 
istic  ice  edge  is  primarily  caused  by  inadequate  tuning  of  the  dynamic  and  thermodynamic 
parameters.  This  possibility  is  supported  by  the  1-year  sensitivity  studies.  The  improvement 
yielded  by  these  simulations  is  illustrated  in  Figure  13,  which  shows  August  thickness  con¬ 
tours,  and  Figure  14,  which  shows  mid-August  compactness  transects.  The  initial  conditions 
for  these  simulations  were  the  31  December  data  from  the  fourth  year  of  the  standard 
simulation.  Consequently,  the  different  ice  edge  characteristics  only  have  about  8  months  to 
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Figure  14.  Mid- August  compactness  transects 
(through  grid  cells  4  and  8)  for  the  standard  and 
two  sensitivity  simulations. 


evolve.  Even  over  this  short  time  there  is  a  significant  improvement  in  both  the  ice  edge  loca¬ 
tion  and  the  amount  of  open  water  formed  in  the  central  Arctic.  Note  also  that  both  sensitiv¬ 
ity  simulations  yield  a  sharper  ice  edge  with  a  rapid  rise  to  about  40%  ice  concentration. 

While  the  increase  in  ice  stress  in  the  “strength”  study  is  probably  excessive,  the  change  in 
parameters  for  the  “growth”  simulation  is  within  the  range  of  uncertainty  of  the  forcing  pa¬ 
rameters.  The  ice  albedo  under  melting  conditions  for  the  high-growth  case,  for  example,  is 
equal  to  Semtner’s  (1976)  adjusted  value  for  multiyear  ice.  While  this  value  may  be  a  bit  high 
for  ice  covered  with  melt  ponds  (see,  e.g.,  Langleben  1971),  summer  radiation  values  are  not 
known  precisely  because  of,  for  example,  uncertainties  in  cloud  cover.  Also,  adding  in  in¬ 
creased  ice  growth  rates  from  snowfall  could  largely  offset  the  oceanic  heat  flux  term  in  the 
simple  thermodynamic  ice  model  used  here. 

Analysis  of  these  sensitivity  results  shows  that  the  reasons  for  the  improvement  of  the  ice 
edge  are  quite  different  for  the  two  sensitivity  simulations.  In  the  case  of  the  thermodynamic 
sensitivity  study,  it  is  mainly  a  case  of  reduced  summer  melting,  which  simply  melts  less  ice. 
In  the  dynamics  case  the  effect  is  primarily  one  of  reducing  offshore  advection,  especially  at 
a  critical  time  in  the  early  spring.  The  advection  effect  is  illustrated  in  Figure  1 5a,  which 
shows  the  difference  between  ice  growth  and  actual  ice  thickness  at  grid  cell  4.  A  major  dif¬ 
ference  between  the  simulations  occurs  in  the  spring  at  approximately  day  120.  At  this  time 
both  the  high-growth  simulations  show  a  sharp  reduction  in  thickness  owing  to  the  dy¬ 
namics,  whereas  the  high-strength  case  does  not.  Examination  of  the  high-strength  case 
shows  that  the  basin  is  practically  frozen  by  the  high  ice  stresses  in  early  spring,  thus  prevent¬ 
ing  such  offshore  motion.  This  lack  of  offshore  motion  in  turn  keeps  the  ice  from  thinning 
early  in  the  melt  season  and,  hence,  reduces  the  amount  of  offshore  open  water.  While  such 
motionless  situations  are  unrealistic  for  this  year  (see  the  next  section),  they  do  occur  under 
certain  onshore  or  quiescent  wind  conditions  as  documented  by  Thorndike  and  Colony 
(1980)  for  the  Beaufort  Sea. 

The  role  of  reduced  summer  melt  for  the  high  growth  case  is  illustrated  in  Figure  15b, 
which  compares  ice  thicknesses  versus  time  at  grid  cell  6.  This  figure  also  demonstrates  the 
importance  of  the  boundary  layer  heat  storage  on  fall  ice  freezeup.  Since  the  open  water  al¬ 
bedo  has  not  been  changed,  this  reduced  summer  melt  is  primarily  caused  by  reduced  ice 
melt.  (A  smaller  portion  of  the  reduction  is  attributable  to  setting  the  oceanic  flux  term  equal 
to  zero.)  On  the  basis  of  this  sensitivity  study  it  seems  likely  that  inclusion  of  a  snow  cover 
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a.  Ice  thickness  and  the  net  ice  growth  at  grid  cell  4. 
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b.  Ice  thickness  and  boundary  layer  heat  stor¬ 
age  (in  equivalent  ice  thickness  units)  at  grid 
cell  6. 


Figure  IS.  Time  series  of  the  difference  between  parameters. 


would  improve  the  simulated  ice  edge.  With  a  snow  cover  it  might  take  several  weeks  to  melt 
the  snow  and  reduce  the  surface  albedo  to  snow-free  values.  In  the  simplified  model  used 
here,  the  surface  albedo  changes  immediately  after  melting  conditions  commence. 

It  is  also  likely  that  some  of  the  excessive  summer  melt  in  the  standard  case  is  caused  by  an 
unrealistically  small  albedo  of  open  water  (0.1).  Observational  estimates  (Langleben  1971) 
are  closer  to  0.2,  a  value  that  is  commensurate  with  albedos  in  the  polar  regions  used  in  most 
climate  models  (see,  e.g.,  Manabe  and  Stouffer  1980). 
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It  is  possible  that  these  sensitivity  results  may  have  relevance  to  ice  edge  forecasting. 
Rogers  (1978),  for  example,  has  found  the  air  temperature,  in  the  form  of  thawing  degree 
days,  is  the  parameter  that  correlates  most  highly  with  the  summertime  ice  margin  off  the 
North  Slope.  Barnett  (1980),  on  the  other  hand,  has  found  that  pressure  changes  in  April 
correlate  well  with  a  combined  severity  index  of  North  Slope  ice  conditions  in  August  and 
September.  The  simulation  done  here  suggests  that,  while  thermodynamics  probably  con¬ 
trols  the  overall  degree  of  melt,  dynamical  effects  early  in  the  spring  may  be  critical  in  thin¬ 
ning  out  the  near-shore  ice.  This  thinning  can  then  produce  a  local  ice  edge  even  under  rela¬ 
tively  “cool”  conditions. 


Ice  thickness  characteristics  off  the  Canadian  Archipelago 

The  ice  thickness  characteristics  off  the  Canadian  Archipelago  also  require  several  years  to 
fully  evolve.  However,  this  evolution  is  one  of  thickness  buildup  rather  than  decay.  Some  of 
the  evolution  characteristics  are  illustrated  in  Figure  16,  which  shows  the  ice  thickness  and 
fraction  of  thick  ice  on  a  transect  between  grid  cells  2  and  8.  Also  shown  in  this  figure  is  the 
total  amount  of  ice  ridged  per  year.  The  observed  data  were  taken  from  submarine  sonar  es¬ 
timates  of  ice  thickness  and  the  thick  ice  fraction  obtained  by  Wadhams  (1981)  using  Octo¬ 
ber  1976  data  from  the  HMS  Sovereign. 

The  main  characteristic  of  the  evolution  is  a  gradual  increase  in  ice  thickness  near  the 
coast,  coupled  with  a  decrease  at  the  Pole.  The  slow  buildup  comes  from  to  new  thin  ice  be¬ 
ing  ridged  each  year,  thus  gradually  increasing  the  thickness.  The  thinning  of  ice  near  the 
Pole,  on  the  other  hand,  appears  to  be  primarily  an  advection  process,  with  thinner  ice  from 
the  Beaufort  Sea  being  brought  in.  The  2-year  lag  is  needed  for  this  thinning  because  the  ice 
in  the  Beaufort  Sea  is  rather  thick  at  First  because  of  the  initial  conditions.  This  decrease  in 
thickness  is  also  accompanied  by  a  more  peaked  and  less  asymmetric  thickness  distribution 
at  the  Pole  on  year  5  than  on  year  2.  This  characteristic  is  illustrated  in  Figure  17,  which 
shows  the  interannual  evolution  of  the  simulated  ice  thickness  distribution  near  the  Pole.  For 
a  comparison  to  observations.  Figure  17  also  shows  observed  ice  drift  distributions  near  the 
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Figure  16.  Evolution  of  the  simulated  ice  thickness 
characteristics  off  the  Canadian  Archipelago  along  a 
transect  through  grid  cells  2  and  8.  The  observed  thick¬ 
ness  results  were  obtained  in  (a)  by  multiplying  ice  draft  data 
from  Wadhams  (1981)  by  U  to  approximately  convert  ice 
drafts  to  thickness  and  in  (b)  by  using  Wadham ’s  areal  per¬ 
centages  for  ice  drafts  deeper  than  5  m. 
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c.  Net  ice  ridged  (in  metres)  over  an  annual  cycle. 

Figure  16  (cont’d).  Evolution  of  the  simulated  ice 
thickness  characteristics  off  the  Canadian  Ar¬ 
chipelago  along  a  transect  through  grid  cells  2  and  8. 
The  observed  thickness  results  were  obtained  in  (a)  by 
multiplying  ice  draft  data  from  Wadhams  (1981)  by  1.1  to 
approximately  convert  ice  drafts  to  thickness  and  in  (b)  by 
using  Wadham 's  areal  percentages  for  ice  drafts  deeper  than 
5  m. 


Pole  obtained  by  Wadhams  (1981)  and  by  LeShack.*  The  LeShack  data  were  obtained  from 
the  February  1960  cruise  of  the  USS  Sargo. 

Figures  16  and  17  show  the  simulated  results  to  have  a  thick  ice  percentage  similar  to  Wad- 
hams’s  data,  but  somewhat  smaller  thickness  values,  especially  near  the  Pole.  Part  of  this 
discrepancy  may  arise  from  insufficient  growth.  In  the  growth  sensitivity  study,  for  example, 
there  is  a  considerable  increase  in  thickness  after  only  1  year.  This  feature  is  shown  in  Figure 
18,  which  compares  sensitivity  transects  off  the  Canadian  Archipelago  with  the  standard 
case.  However,  there  may  also  be  significant  amounts  of  interannual  variability.  The 
LeShack  data  at  the  Pole,  for  example,  are  in  better  agreement  with  the  simulated  values 
than  Wadhams’s  data. 

An  interesting  feature  shown  by  the  net  ridging  is  a  slow  evolution  of  a  zone  of  high  ridg¬ 
ing  several  hundred  kilometres  off  the  coast.  This  feature  does  not  occur  the  first  few  years, 
which  indicates  that  it  is  probably  related  to  the  spatial  variations  in  ice  strength  that  are 
gradually  built  up.  While  this  ridging  peak  is  not  manifested  in  the  simulated  thick  ice  per- 

•Personal  communication  with  L.A.  LeShack,  LeShack  Associates,  Silver  Spring,  Maryland,  1979. 


26 


0  32 


Initial 

Conditions 


§p 

I  Stir* I 


Yeor  2 ,  day  129 

b. 


Thickness  <m) 


Year  5,  day  129 

C. 


GRID  CELL  8 


llrflP 


Thickness  (m) 


i 


Sovereign 


#  •• 


('4® 


■mt&m 


GRID  CELL  8 


Thickness  (m) 


m  m 

,•  -  .V;’  1 

**fki  ’  L~1— ? 


fW&m 


Ice  Droft  (m) 


Figure  17.  Interannual  evolution  of  the  simulated  ice  thickness  distribution  at  the  Pole  (grid 
cell  8).  The  observed  ice  draft  distributions  were  obtained  by  Wadhams  (1981)  using  October  1976 
HMS  Sovereign  data,  and  by  LeShack  (persona!  communication)  using  February  1960  USS  Sargo  data. 


centage  (which  falls  off  smoothly)  there  is  an  indication  of  a  discontinuous  fall-off  at  about 
500  km  in  Wadhams’s  observed  thick  ice  percentages.  Also,  and  perhaps  more  relevant, 
basin-wide  ice  roughness  contours  compiled  by  LeShack*  indicate  such  an  offshore  rough¬ 
ness  peak  in  this  region. 

The  relative  roles  of  ridging  and  growth  in  maintaining  the  thick  ice  near  the  coast  are 
demonstrated  in  Figure  19.  This  Figure  shows  the  net  ice  transfer  by  ridging  and  vertical 
growth  at  grid  cell  2  over  an  annual  cycle  together  with  the  thickness  distribution  on  day  129. 


'Personal  communication  with  L.A.  LeShack,  LeShack  Associates,  Silver  Spring,  Maryland,  1979. 
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Figure  18.  Ice  thickness  transects  through  grid  cells  2 
and  8  for  the  standard  and  the  two  sensitivity  simula¬ 
tions. 
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Figure  19.  Ice  thickness  distribution  (a)  at  grid  cell  2,  and 
net  annual  ice  transfer  probability  density  (b)  owing  to 
growth  and  ridging  at  grid  cell  2. 
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Over  the  year,  ridging  transfers  thin  ice  up  to  about  3  m  thick  to  thicker  categories.  The 
growth,  on  the  other  hand,  causes  a  net  loss  of  thick  ice.  In  seasonal  equilibrium,  imbalances 
between  these  two  terms  are  accounted  for  by  advection  and  deformation  (and  to  a  much 
lesser  degree  lateral  melt).  For  comparison  it  is  worth  noting  that  plots  similar  to  these  yield 
ice  up  to  only  about  2  m  being  ridged  near  the  North  Pole  and  ice  up  to  only  about  1  m  being 
ridged  in  the  southern  Beaufort  Sea. 

Comparison  of  observed  and  simulated  ice  drift 

For  a  direct  comparison  to  observation,  the  simulated  drift  rates  of  three  drifting  ice  sta¬ 
tions  were  determined  by  interpolation.  Table  1  compares  some  of  the  simulated  deforma¬ 
tion  rates  and  net  drift  to  the  observed  values.  The  net  drift  is  the  average  drift  from  day  140 
(1962)  to  day  109  (1963)  of  the  three  drifting  ice  stations.  The  x  and  y  strain  rates  are  based 
on  11 -day  averaged  velocity  fields  at  the  ice  station  locations.  Except  for  the  high-strength 
case,  which  yields  small  deformation  fluctuations,  these  comparisons  show  relatively  similar 
statistical  behavior  for  the  simulated  deformation  rates  but  significant  differences  in  the 
long-term  drift.  Specifically,  the  simulated  net  drift  for  the  high  strength  case  is  almost  half 
that  of  the  standard  case.  It  is  also  notable  that  in  the  standard  simulation  the  fifth-year  net 
drift  is  about  30Vo  larger  than  the  first-year  because  of  a  gradual  reduction  in  ice  strength 

Much  of  the  reduction  in  net  drift  in  the  high-strength  case  is  caused  by  an  effective  stop¬ 
page  of  almost  all  motion  in  April  and  May.  This  feature  is  apparent  from  Figure  20,  which 


Figure  20.  Simulated  and  observed  drift  rates  of  ice  station  Artis  f May  1962-April  1963). 


Table  1.  Observed  and  simulated  ice  station  drift. 


5th  year 
standard 
simulation 

1st  year 
standard 
simulation 

High- 

strength 

simulation 

High- 

growth 

simulation 

Observed 

Net  ice  station  drift 

702  km 

567  km 

339  km 

682  km 

562  km 

Net  drift  angle 

26  °W 

27  °W 

13°W 

26  °W 

6°W 

Standard  deviations 
of  strain  rate: 

*xx 

0.0033  day'1 

0.0038  day’1 

0.0028  day' 

0.0035  day-1 

0.0028  day' 

fyy 

0.0024  day-' 

0.0024  day'1 

0.0021  day-' 

0.0025  day-' 

0.0033  day"1 

Correlation  coeffi¬ 
cients  between  simu¬ 
lated  and  observed 
strain  rates 

0.33 

0.38 

0.31 

0.36 

(yy 

0.65 

0.65 

0.56 

0.70 

- 

Figure  21.  Outflow  time  series  for  the  standard  and  high-strength  simulations.  (Units  are 
metres  of  ice  averaged  over  the  basin.  For  comparison  1.0  m/yr  is  equivalent  to  0.225  Sv.) 


compares  the  simulated  and  observed  drift  rates  for  ice  station  Arlis.  The  ice  velocity  also 
stops  in  the  outflow  region  between  Greenland  and  Spitsbergen,  as  demonstrated  by  the  ice 
thickness  outflow  time  series  in  Figure  21.  This  “freezing  up”  of  the  basin  takes  place  pri¬ 
marily  in  response  to  art  increase  in  the  ice  strength  resulting  from  removal  of  thin  ice  by 
growth  and  deformation.  A  secondary  effect  may  be  changes  in  the  mean  wind  patterns  as 
the  high  pressure  system  over  the  Beaufort  Gyre  weakens  and  the  ice  circulation  is  more 
dominated  by  a  low  pressure  system  approximately  centered  over  Severnaya  Zemlya  (Soviet 
islands  at  upper-right-hand  portion  of  grid).  This  low  pressure  pattern  is  most  pronounced  in 
March  when  the  stoppage  begins.  There  is  also  some  reduction  in  mean  wind  speeds  in  April 
and  May,  but  this  effect  appears  secondary,  as  illustrated,  for  example,  by  the  substantial 
simulated  drift  rates  of  the  observed  and  the  standard  cases  in  April  and  May  (Fig.  20). 

Temporal  variations  in  ice  strength  for  both  the  standard  and  high-strength  simulations 
are  shown  in  Figure  22.  The  seasonal  variation  in  strength  occurs  in  response  to  the  thermo¬ 
dynamic  forcing  which  in  summer  melts  ice  and  in  winter  removes  thin  ice  by  freezing.  This 
seasonal  variation  is  consistent  with  seasonally  varying  best-fit  ice  strengths  obtained  by  Hib- 
ler  and  Tucker  (1977)  and  with  analysis  of  ice  drift  in  summer  by  McPhee  (1980b).  Note  that 
the  main  reduction  in  strength  takes  place  about  at  day  160  when  ice  begins  melting  in  the 
central  basin  (see  Fig.  6).  However,  in  the  standard  case  some  seasonal  weakening  begins  ear¬ 
lier.  This  probably  happens  because  open  water  created  by  deformation  is  not  removed  by 
freezing  after  about  day  120.  Superimposed  on  this  seasonal  variation  are  fluctuations  that 
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Figure  22.  Ice  strength  time  series  for  the  standard 
and  high-strength  simulations. 
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are  especially  pronounced  in  the  near-shore  region.  These  large  fluctuations  arise  because  tne 
ice  strength  is  determined  only  by  the  thinnest  15*7o  of  the  ice  cover  (see  Appendix  A). 

It  is  notable  that  the  standard  case  strengths  are  about  an  order  of  magnitude  higher  than 
those  obtained  by  Thorndike  et  al.  (1975).  While  part  of  this  can  be  attributed  to  including 
frictional  losses,  the  main  increase  is  caused  by  the  way  the  redistributor  used  here  transfers 
thin  ice  into  thicker  categories  under  deformation.  It  is  also  notable  that  the  strength  range 
spanned  by  the  standard  and  high-strength  simulations  is  similar  to  the  range  deduced  from 
near-shore  ice  dynamics  studies  by  Pritchard  (1978).  In  particular,  using  a  spatially  constant 
strength  in  a  localized  near-shore  region,  Pritchard  found  a  best-fit  compressive  strength  of 
*4 x  104  N  m'1.  When  the  compressive  strengths  were  an  order  of  magnitude  higher  than 
this,  the  ice  velocity  field  was  almost  entirely  determined  by  the  boundary  motion.  The  simu¬ 
lations  conducted  here  indicate  that  such  a  dominance  by  boundary  motion  can  also  occur  in 
a  fully  coupled,  basin-wide  dynamic-thermodynamic  model  when  frictional  losses  in  ridging 
are  sufficiently  increased.  This  result  may  have  relevance  to  paleoclimate  arguments  con¬ 
cerning  stoppage  of  the  ice  flow  from  the  Arctic.  In  particular,  these  results  suggest  that 
higher  growth  rates  coupled  with  a  modest  increase  in  frictional  losses  might  cause  the  basin 
to  effectively  freeze  up  in  winter. 


Mtu  balance  characteristics 

A  dominant  feature  of  the  variable  thickness  model  is  an  increased  seasonal  variation  in 
the  ice  growth  and  decay.  This  feature  is  illustrated  in  Figure  23,  which  shows  the  basin-aver¬ 
aged  thickness  time  series  of  the  standard  and  thermodynamics-only  simulations.  While 
some  of  the  increased  growth  is  ascribable  to  a  large  amount  of  open  water  in  the  fall,  much 
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Fcr  computation,  eq  B1  may  be  rewritten  in  the  form: 

Al+A1(To)  +  (K/H)(T„-T0)  =  0  (B2) 

where  A,  is  a  constant  and  A2(T0)  is  a  nonlinear  function  of  T0.  Using  a  Newton-Raphson 
procedure,  one  solves  eq  B2  for  T0.  If  T0  is  greater  than  273.16  K  it  is  set  equ-xl  to  273.16  K.  In 
both  cases  the  growth  rate  is  then  calculated  from 

f,(H)  =  -  [/I,  +  A 2(To)  +  F0]  /  Q,  (B3) 

where  Q,  is  the  volumetric  heat  of  fusion  of  ice,  set  at  3.02  x  10'  J  m~\  and  F0  is  a  constant 
oceanic  heat  flux  into  the  mixed  layer  from  the  deep  ocean.  In  the  case  of  open  water,  eq  B3 
is  also  used  to  estimate  the  growth  rate,  with  the  modification  that  the  albedo  is  set  equal  to 
0.1,  To  is  set  equal  to  Tw,  and  the  latent  heat  transfer  coefficient  D,  is  slightly  different.  In  any 
given  grid  cell  the  water  temperature  is  calculated  from  the  mixed  layer  equations  discussed 
earlier  (eq  6  and  7),  and  in  the  presence  of  the  ice  is  always  equal  to  271.2  K.  Note  that  while 
snow  cover  is  not  considered  explicitly  in  this  model  it  is  easy  to  determine  the  growth  rate  of 
ice  with,  say,  a  snow  cover  of  depth  hs  on  top.  For  this  purpose  (following  Semtner  [1976]), 
K  in  eq  Bl  and  B2  is  simply  replaced  by 

K, 

(h,/H)  +  (K,/K) 

where  K,  is  the  snow  conductivity. 

To  examine  the  behavior  of  this  heat  budget  code  in  more  detail,  the  individual  compo¬ 
nents  were  compiled  every  16  days  at  selected  locations.  Table  Bl  lists  the  components  of  the 
heat  balance  at  grid  cell  8  for  winter  growth  conditions  (18  February),  summer  melt  condi¬ 
tions  (26  June),  and  for  a  growth-melt  transition  period  (10  June).  The  atmospheric  input 
values  for  these  times  and  locations  are  shown  in  Table  B2.  For  comparison,  also  shown  in 
these  tables  are  central  Arctic  heat  balance  estimates  obtained  by  Maykut  (1978)  together 
with  his  forcing  components. 

In  making  this  comparison  it  should  be  noted  that  in  addition  to  different  forcing  fields 
there  are  differences  between  Maykut’s  (1978)  sea  ice  thermodynamic  formulation  and  that 
used  here.  In  particular,  for  thin  ice  up  to  a  metre  thick,  Maykut  allows  the  conductivity  to 
be  temperature-dependent  and  the  surface  albedo  to  decrease  with  thickness.  He  does,  how¬ 
ever,  assume  a  steady  state  equilibrium  heat  budget.  For  thick  ice  (3  m),  Maykut  uses  the 


Table  B2.  Central  Arctic  (grid  cell  8)  atmospheric  tent' 
peratures,  dew  points  and  wind  speeds.* 


18  February 

10  June 

16  June 

Air  temperature  (K) 

241 

269 

272 

(242) 

(269) 

Wind  speed  (m  s'1) 

6.17 

6.18 

1.47 

(5) 

(5) 

Dew  point  temperature  (K) 

236 

268 

271 

'Central  Arctic  values  used  by  Maykut  are  given  in  parentheses  and 
represent  conditions  on  the  first  day  of  the  the  month.  Maykut  as¬ 
sumed  a  relative  humidity  of  0.90  in  February  and  0.96  in  June. 
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APPENDIX  B:  HEAT  BUDGET  CODE 


To  solve  for  the  ice  growth  rate  a  surface  heat  budget  computation  similar  to  that  of  Man- 
abe  et  al.  (1979)  and  Parkinson  and  Washington  (1979)  is  used.  This  heat  budget  code  incor¬ 
porates  a  time-independent  thermodynamic  sea  ice  model  similar  to  the  simplest  of  a  hier¬ 
archy  of  models  developed  by  Semtner  (1976).  Semtner’s  simplest  time-independent  model 
utilizes  a  constant  conductivity,  together  with  a  linear  temperature  profile.  In  addition,  a 
seasonal  snow  cover  is  allowed  to  build  up  and  decay.  Both  the  conductivity  and  the  snow- 
free  surface  albedo  were  adjusted  by  Semtner  (1976)  to  compensate  for  not  including  internal 
melting  in  summer.  With  these  adjustments,  Semtner  was  able  to  produce  mean  annual 
thicknesses  that  agreed  very  well  with  a  more  complete  time-dependent  model  developed  by 
Maykut  and  Untersteiner  (1971).  In  the  computations  done  here,  Semtner’s  (1976)  time-inde¬ 
pendent  approach  is  used.  However,  following  Manabe  et  al.  (1979)  the  effects  of  snow 
cover  are  approximated  only  implicitly  by  changing  the  surface  albedo  of  the  ice  to  snow  val¬ 
ues  for  below-freezing  surface  temperatures. 

The  basic  surface  heat  balance  equation  in  the  presence  of  an  ice  cover  is  (where  fluxes  into 
the  surface  are  considered  positive) 

(1  - a)Fs  +  FL  +  D,  \Uc\  (r.-r„)  +  A  |{/c|  [q.(Z)-q,(To)] 

-  D, Ti  +  (K/H)(T. -To)  =  0  (B 1 ) 

where  a  =  surface  albedo 

To  =  surface  temperature  of  the  ice 
T.  =  air  temperature 
K  =  ice  conductivity 

H  -  ice  thickness  (which  for  computation  is  set  equal  to  0.05  m  for  very  thin  ice) 

Tw  =  water  temperature 

Uc  =  geostrophic  wind 

q.  -  specific  humidity  of  the  air 

q,  =  specific  humidity  at  the  ice  surface 

F,  =  incoming  shortwave  radiation 

Fi  -  incoming  longwave  radiation. 

The  constants  D ,  and  D%  are  bulk  sensible  and  latent  hr  at  transfer  coefficients  and  D,  is  the 
Stefan-Boltzmann  constant  times  the  surface  emissivity.  Numerical  values  identical  to  those 
used  by  Parkinson  and  Washington  (1979)  are  used  for  these  constants.  (Specifically,  D,  = 
2.28  J  m'J  K'1,  D,  =  5.50 x  10‘*  W  nrJ  K'*,  A  =  5.69  x  10’  J  m'1  over  water  and  6.45  x  10’  J 
nr’  over  ice.)  The  specific  humidities  at  the  ice  (and  water)  surfaces  are  calculated  in  the 
same  manner  as  described  by  Parkinson  and  Washington  (1979).  Specific  humidities  and  air 
temperatures  at  a  nominal  10-m  reference  level  above  the  ice-water  surface  were  obtained 
from  climatological  data  compiled  by  Crutcher  and  Meserve  (1970).  For  the  ice  albedo  0.616 
was  used  when  the  surface  temperature  was  equal  to  273.16  K.  This  value  is  slightly  smaller 
than  Semtner’s  adjusted  value  of  0.66  but  larger  than  Parkinson  and  Washington’s  adjusted 
value  of  0.53  for  snow-free  ice.  For  ice  with  a  surface  temperature  below  freezing,  a  =  0.75 
was  used,  which  is  identical  to  Parkinson  and  Washington’s  “snow”  albedo.  For  the  con¬ 
ductivity  Semtner’s  adjusted  value  of  2.1656  W  m'1  K'1  was  used. 
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Figure  Al.  Ice  strength  versus  thickness  for  the 
redistributor  proposed  here.  Strengths  for  different 
values  of  the  parameter  H*  (in  metres )  are  shown. 


highly  dependent  on  2 H*.  Consequently,  uncertainties  in  H *  are  not  critical  to  the  simula¬ 
tions.  It  is  notable  that  strengths  are  substantially  higher  than  those  obtained  by  Thorndike 
et  al.  In  particular,  for  25-cm-thick  ice,  the  Thorndike  et  al.  (1975)  redistributor  would  yield 
a  strength  of  about  5  x  10J  N  m"1  whereas  this  redistributor  yields  =2.5  x  103  N  m'1. 

It  should  be  noted  that  there  are  some  similarities  between  this  redistributor  and  that  pro¬ 
posed  by  Bugden  (1979)  and  that  used  in  early  developments  of  the  ice  thickness  distribution 
theory  (e.g.  Thorndike  and  Maykut  1973).  Bugden’s  contribution  was  to  show  that  such  a 
constant  redistributor  could  at  least  partially  be  deduced  from  ridge  statistical  relations  ob¬ 
tained  by  Hibler  et  al.  (1972)  and  Mock  et  al.  (1972).  The  original  feature  of  the  redistribu¬ 
tion  proposed  here  is  that  the  square  root  scaling  of  the  maximum  thickness  cutoff  be  based 
on  geometric  arguments. 


where  c>  is  a  constant,  say  0.  IS.  The  basic  idea  here  is  that  under  ridging  conditions  the  clos¬ 
ure  of  both  thin  ice  and  thick  ice  occurs;  however,  the  thinnest  ice  is  allowed  to  be  removed 
to  a  greater  extent.  By  fixing  c,  to  be  0.15,  only  the  thinnest  15%  of  the  ice  is  deformed.  Note 
that  if  there  is,  say,  25%  open  water,  only  open  water  will  be  removed  and  there  will  be  no 
ridging.  In  the  author’s  judgment  this  choice  of  P(h)  is  quite  reasonable  and  is  adopted  here 
with  c,  =  0.15. 

The  other  unknown  is  y(h,,  hi),  which  specifies  how  ice  is  transferred  from  one  category  to 
another  by  ridging.  Thorndike  et  al.  (1975)  have  suggested  that 

>(*.,*»)  =  5<A,-Ar/i.Kl/Ar)  (A9) 
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where  k  is  a  constant.  Equation  A8  effectively  states  that  ice  is  transferred  into  a  fixed  multi¬ 
ple  of  its  own  thickness.  This  form  is  computationally  simple  and  is  useful  for  examining  the 
behavior  of  the  ice  thickness  distribution  model.  However,  the  redistribution  of  eq  A8  ig¬ 
nores  two  important  physical  features  of  the  ridging.  One  feature  is  that  under  deformation, 
ridging  transfers  ice  to  a  variety  of  thickness  categories.  This  is,  for  example,  evident  from 
field  observations  of  ridges  (see,  e.g..  Weeks  et  al.  1971,  and  Kovacs  1972),  which  show  them 
to  be  roughly  triangular  in  shape.  These  observations  indicate  that  under  deformation,  leads 
containing  thin  ice  of  a  given  thickness  typically  form  ridges  having  a  triangular  cross  sec¬ 
tion.  To  satisfy  eq  A9  such  ridges  would  have  to  have  vertical  sides. 

The  second  physical  feature  is  that  typical  ridge  heights  divided  by  the  thickness  of  ice  be¬ 
ing  ridged  appear  to  decrease  with  increasing  thickness  (see,  e.g..  Tucker  and  Govoni  1981). 
This  particular  scaling  is  an  important  feature  of  the  Parmerter  and  Coon  (1972)  mechanistic 
ridge  model.  In  particular,  Parmerter  and  Coon  (1972)  found  that  ridges  simulated  by  using 
their  mechanistic  model  tended  to  have  a  maximum  limiting  height.  This  limiting  height, 
however,  tended  to  scale  approximately  as  the  square  root  of  the  thickness  of  ice  being  ridged. 
Such  a  square  root  dependence  is  also  consistent  with  one’s  intuition  concerning  ridging. 
Consider,  for  example,  two  equal  width  leads  undergoing  ridging.  With  the  assumption  that 
all  the  thin  ice  is  deformed  into  piles  of  relatively  small  blocks  with  similar  slope  angles,  the 
two  ridges  formed  would  have  heights  scaling  with  the  square  root  of  the  thickness  of  the  ice 
being  ridged.  If  the  distribution  of  lead  widths  is  independent  of  the  ice  thickness  in  leads, 
such  an  intuitive  argument  should  have  application. 

A  number  of  redistributors  could  be  constructed  that  approximately  satisfy  these  con¬ 
straints.  However,  the  simplest  case  is  to  take  y  to  be  a  constant  up  to  some  cutoff  thickness: 


y(h,,  hi)  =  Ci  for  2 h,  s  hi  s  b(/r,)A, 


where 


2(/f*/h,)“ 


with  H*  a  constant.  Examination  of  Parmerter  and  Coon’s  (1972)  mechanistic  studies  sug¬ 
gests  that  H*  -  100  m  is  a  reasonable  value.  For  1-m-thick  ice  this  yields  a  maximum  ridged 
ice  thickness  of  20  m.  What  type  of  strengths  this  yields  is  illustrated  in  Figure  Al,  which 
shows  typical  ice  strengths  (obtained  numerically)  versus  thickness  for  this  redistributor.  For 
comparison,  strengths  for  several  different  values  of  H*  are  shown.  (Since  ice  thicker  than 
18.3  m  is  not  allowed  numerically,  the  thick  ice  strengths  in  Figure  Al  will  be  somewhat 
smaller  than  analytical  values  obtained  in  the  Analytic  Examination  of  the  Ridge  Redistribu¬ 
tion  Process  section.)  An  item  of  particular  interest  for  this  paper  is  that  the  strengths  are  not 
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(A5) 


The  first  term  in  eq  A4  specifies  the  amount  of  open  water  created,  while  the  second  term  de¬ 
scribes  the  transfer  of  thin  ice  to  thicker  categories  by  ridging.  In  this  formalism  y(hlt  h2)dh2 
can  be  thought  of  as  the  area  of  ice  put  into  the  thickness  interval  [h2,  h2  +  dh2]  when  a  unit 
area  of  ice  of  thickness  h,  is  used  up.  Basically,  the  integral  over  y  specifies  where  the  ice  is 
transferred  to  during  ridging.  The  function  P(h),  on  the  other  hand,  is  some  probability 
function  specifying  which  categories  are  being  destroyed  by  ridging.  Written  in  this  form  7 
automatically  satisfies  eq  Al.  Satisfaction  of  eq  A2  and  A3  leads  to  the  additional  con¬ 
straints,  respectively. 


f  y(h h)hdh 


0 


=  h ' 


(A6) 


OO 

C  J  WM'dh  =  P*.  (A7) 

0 

Equation  A6  requires  that  transfer  of  ice  from  the  category  does  not  change  its  mass,  and  eq 
A7  serves  as  a  definition  of  the  ice  strength  P*. 

To  get  some  feeling  for  the  amount  of  open  water  created  with  this  formulation,  it  is  use¬ 
ful  to  explicitly  calculate  the  coefficients  in  eq  A4.  The  dynamical  code  used  here  employs  an 
elliptical  yield  curve  (see  Hibler  1979).  Using  this  particular  rheology  (and  noting  that  the 
symbol  P  in  eq  4  of  Hibler  [1979]  is  equivalent  to  the  symbol  P*  used  here)  one  can  easily 
show  that: 


KeJ/P*  =  0.5(A-e„) 


where 

A  =  [(ef, +€10(1  +  1/0  +  4*-J  €12  +  2e, ,  -  l/e1)]" . 

Note  that  A  is  dependent  only  on  the  two  invariants  of  the  strain  rate  tensor:  i,  =  (€,,+€22) 
and  €«  =  [(€22  -  €,  O1  +  4e?a].  In  this  expression,  e  is  the  ratio  of  the  principal  axes  of  the  ellip¬ 
tical  yield  curve  (set  equal  to  2)  and  is  a  measure  of  the  relative  shear  strength  to  compressive 
strength.  Larger  e  values  yield  less  shear  strength.  With  this  rheology  no  open  water  will  be 
formed  under  pure  convergence  (€11  —  €h\  €12  —  0;  €  1 1  <  0).  Conversely,  under  pure  diverg¬ 
ence  (f,  1  =  €22;  €12  =  0;  €n  >  0)  there  will  be  no  ridging.  Also,  as  e  —  00  the  amount  of  open 
water  forming  under  pure  shear  (e,,  =  €22  =  0;  €,2  #  0)  will  approach  zero.  In  the  Thorndike 
et  al.  (1975)  study  the  special  case  of  e  =  1  was  used. 

The  two  main  unknowns  in  the  redistribution  theory  are  P(/>)  and  7 (h,,  h2).  Thorndike  et 
al.  (1975)  have  suggested  that  P{'  )  be  given  by 

P[h)  =  Max  [(-  «(/r)dh/c^,oJ  (A8) 
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APPENDIX  A:  MECHANICAL  REDISTRIBUTOR 


Consider  a  two-dimensional  continuum  with  ice  that  is  described  by  a  plastic  rheology 
such  that  the  stress  state  a0(e0t  P*)  is  a  function  of  the  strain  rate  e0  and  a  strength  P*.  Then, 
following  Thorndike  et  al.  (1975)  and  Rothrock  (1975),  the  mechanical  redistribution  func¬ 
tion  (see  eq  4)  must  satisfy  the  following  constraints: 


w 

J*  ^  dh  =  V  •  n 


j  h*dh  =  0 


QD 

C  J  W(A)d!ft  = 


where  u  is  the  ice  velocity  field  and  C  is  a  constant.  In  consideration  of  only  the  work  done 
on  the  ice  by  gravitational  buoyancy  forces,  C  would  be  given  by 


c»  =  'A  e,|[(e„ -&)/<?.$} 


with  q,  the  density  of  ice,  e»  the  density  of  water,  and  ^the  acceleration  of  gravity.  Equation 
Al  follows  from  the  constraint  that  i  renormalize  the  g  distribution  to  unity  because  of 
changes  in  area.  Equation  A2  follows  from  conservation  of  mass  and  basically  states  that  \p 
does  not  create  or  destroy  ice  but  merely  changes  its  distribution.  (An  additional  assumption 
in  eq  A2  is  that  the  ice  mass  is  related  in  a  fixed  manner  to  the  thickness.)  The  third  con¬ 
straint  (Rothrock  1975)  specifies  that  work  done  in  building  ridges  is  equal  to  the  deforma- 
tional  work.  An  important  feature  of  this  constraint  is  that,  for  an  arbitrary  plastic  yield 
curve,  it  forces  some  ridging  to  take  place  even  though  there  may  be  no  net  convergence.  In 
its  present  form  (with  C  =  Ct)  eq  A3  ignores  frictional  losses  occurring  in  ridging  and  con¬ 
siders  only  the  potential  energy  change  from  deformed  ice  being  lifted  against  gravity  and  be¬ 
ing  forced  down  against  buoyancy.  However,  Rothrock  (1975)  has  estimated  these  frictional 
losses  to  be  the  same  order  of  magnitude  as  the  potential  energy  changes.  Consequently,  as  a 
crude  approximation  for  frictional  losses,  eq  A3  is  retained  here  with  a  constant  C  equal  to 
2C*. 

As  in  Thorndike  et  al.  (1975),  a  redistribution  that  satisfies  these  constraints  is  (using  the 
convention  that  repeated  subscripts  are  summed  over) 
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The  interannual  evolution  of  the  summer  ice  edge  is,  on  the  other  hand,  dominated  by 
thermodynamic  growth  and  decay  rates  rather  than  ridge  buildup.  Specifically,  lower  annual 
thermodynamic  net  growth  coupled  with  offshore  ice  advection  gradually  thins  the  ice  off 
the  North  Slope  and  Siberian  coasts.  This  thin  ice  presence,  coupled  with  high  summer  melt¬ 
ing  rates,  creates  an  excessive  ice  edge  in  summer.  Sensitivity  analyses  show  that  the  magni¬ 
tude  of  this  summer  ice  edge  is  critically  dependent  on  the  albedo  used  for  melting  sea  ice.  By 
modifying  this  albedo  by  only  about  10*70,  much  more  realistic  ice  edge  values  may  be  ob¬ 
tained.  Similarly,  small  changes  in  the  open  water  albedo  should  lead  to  substantially  im¬ 
proved  summer  ice  margins.  These  ice  edge  results  emphasize  the  need  for  a  more  detailed 
analysis  of  the  thermodynamic  sensitivity  of  such  variable  thickness  sea  ice  models.  One  ap¬ 
proach  is  to  do  some  process  sensitivity  studies  using  observed  ice  motion  estimates  as  driv¬ 
ing  fields.  Another  approach  would  be  to  compare  observed  and  simulated  ice  edge  results 
over  several  sequential  years. 

The  heat  exchange  characteristics  simulated  here  are  substantially  affected  by  sea  ice  ridg¬ 
ing  and  deformation.  On  a  basin-averaged  scale,  ridging  transfers  thin  ice  to  thicker  cate¬ 
gories,  thus  making  room  for  more  thin  ice  to  form.  The  thicker  ice  thus  formed  grows  slow¬ 
ly  in  winter,  but  ablates  rapidly  in  summer.  Spatial  imbalances  between  ridging  and  open 
water  creation  cause  substantial  spatial  variations  in  the  air-sea  heat  exchange.  Regions  with 
more  ridging  tend  to  have  a  net  heat  gain  from  the  atmosphere.  Regions  of  large  open  water 
creation,  on  the  other  hand,  have  more  growth  and  thus  have  a  net  loss  of  heat  to  the  atmos¬ 
phere.  These  variations  are  in  contrast  to  thermodynamic  simulations,  where  in  equilibrium 
the  net  annual  ice  growth  is  zero  everywhere. 

In  coupled  dynamic-thermodynamic  sea  ice  models,  the  way  in  which  thin  ice  is  statistical¬ 
ly  redistributed  into  thicker  ice  can  significantly  affect  the  geophysical  stresses  in  pack  ice. 
The  agreement  with  ridge  morphological  data  and  the  realistic  simulated  ice  thickness  and 
ridge  buildup  support  the  redistribution  function  proposed  here.  However,  a  major  un¬ 
known  factor  is  the  amount  of  frictional  energy  lost  during  ridging.  The  strengths  simulated 
here  were  realistic,  although  a  bit  small.  But,  sensitivity  analysis  suggests  that  only  a  modest 
increase  in  frictional  losses  during  ridging  would  be  adequate  for  more  realistic  ice  velocities 
and  strengths.  Since  frictional  losses  in  ridging  are  not  precisely  known,  such  increases  are 
physically  justifiable.  Sensitivity  studies  with  more  detailed  ice  velocity  verification  fields 
should  help  determine  necessary  strength  increases  more  precisely.  Polar  drifting  buoy  data 
taken  should  prove  helpful  in  such  studies. 
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Figure  25.  Annual  net  growth  contours  (at  0.6-m  intervals) 
for  the  standard  simulation. 


thin  ice  to  occur  by  divergence  along  one  coast  while  converging  ice  into  thicker  categories  at 
another  coast.  The  thicker  ice  ablates  faster  in  summer  than  it  freezes  in  winter,  while  the  re¬ 
verse  is  true  of  thin  ice.  As  a  consequence,  over  an  annual  cycle  (neglecting  oceanic  effects) 
areas  with  high  ridging  will  have  a  net  heat  gain  from  the  atmosphere  and  divergent  thin  ice 
regions  will  have  a  net  loss. 


CONCLUDING  REMARKS 

This  report  describes  a  variable  thickness  dynamic-thermodynamic  sea  ice  model  and  ex¬ 
amines  the  seasonal  equilibrium  characteristics  of  this  model.  To  develop  this  model,  the 
Thorndike  et  al.  (1975)  ice  thickness  distribution  framework  was  extended  to  include  an  oce¬ 
anic  mixed  layer  at  a  fixed  depth  and  lateral  melting  effects.  In  addition,  a  mechanical  redis¬ 
tribution  process  that  is  consistent  with  hypothesized  and  observed  physics  of  the  ridging 
process  has  been  proposed.  By  combining  this  framework  with  a  viscous-plastic  ice  dynamics 
model  that  was  developed  previously  (Hibler  1979)  and  a  spatially  varying  surface  heat 
budget,  seasonal  equilibrium  simulations  may  be  performed. 

A  dominant  feature  of  the  Arctic  simulations  discussed  here  is  the  time  needed  for  the 
thickness  characteristics  to  fully  evolve.  The  main  characteristic  is  a  buildup  of  ice  along  the 
Canadian  Archipelago  in  conjunction  with  a  thinning  of  ice  off  the  Alaskan  and  Siberian 
coasts.  These  geographical  chi»“  -,  however,  take  several  years  to  fully  develop.  Off  the 
Archipelago,  the  buildup  is  based  upon  an  accumulation  of  ridged  ice  formed  over  several 
years.  During  each  year  the  ice  strength  builds  up  as  thin  ice  is  removed  by  ndging  and 
growth.  This  creates  a  balance  between  the  internal  and  external  stresses.  Since  the  thick  ice 
does  not  directly  affect  the  ice  stresses,  its  equilibrium  balance  requires  a  longer  time  scale. 
Once  equilibrium  is  reached,  summer  ablation  of  the  thick  ice  largely  offsets  the  new  thick 
ice  formed  by  ridging. 


,  \  .V 


.•-Y-v.-.  V\N’ 


'  v.o, 


>  "  .  - 
"  •*  "  -v*  A." 


34 


«*  V 


Figure  24.  Basin-averaged  net  ice  growth  (a)  and  net  ice  trans¬ 
fer  by  ridging  (b)  for  three  categories  of  ice  in  the  standard 
simulation. 

also  be  that  not  enough  open  water  is  being  created  owing  to  inadequate  temporal  and  spatial 
resolution  of  the  wind  fields.  For  example,  while  the  basin-averaged  thin  ice  growth  (1  m)  ap¬ 
proximately  equals  thicker  ice  growth  in  January,  the  ice  growth  in  the  central  basin  in  con¬ 
siderably  less.  From  day  1  to  day  60,  thin  ice  growth  at  grid  cells  7  and  8  accounted  for  only 
about  20%  of  the  ice  growth.  These  values  are  substantially  less  than  estimates  by  Maykut 
(1978),  which  yield  approximately  equivalent  growth  by  thin  and  thick  ice  in  the  central  basin 
in  winter. 

In  addition  to  increasing  the  total  air-sea  heat  exchange  by  ridging  and  open  water  crea¬ 
tion,  the  dynamics  also  causes  a  spatial  imbalance  in  the  heat  exchange.  Because  of  this  im¬ 
balance,  the  net  growth  over  an  annual  cycle  will  not  be  zero  everywhere,  as  in  a  thermo- 
dynamics-only  simulation,  but  will  vary  spatially.  These  variations  can  be  very  significant  as 
demonstrated  in  Figure  25.  This  imbalance  can  occur  by  freezing  ice  at  one  location  and  then 
transferring  it  to  another  where  it  is  melted.  However,  in  practice  the  process  is  more  com¬ 
plex  and  it  relies  upon  ridging  to  a  large  extent.  The  ice  dynamics,  for  example,  can  cause 
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Figure  23.  Basin-averaged  ice  thickness  time  series  for  the  standard  and  ther¬ 
modynamics-only  simulations. 
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is  from  thin  ice  created  by  deformation.  The  increased  decay  rates,  on  the  other  hand,  arise 
because  of  higher  heat  absorption  by  open  water  and  by  loss  of  ice  from  the  basin  by  out¬ 
flow.  Some  of  these  characteristics  are  illustrated  in  Figure  24,  where  the  basin-averaged  net 
growth  and  net  ice  transfer  by  ridging  are  plotted  for  three  broad  ice  categories.  On  a  basin- 
averaged  basis,  open  water  creation  is  balanced  by  ridging,  which  transfers  thin  ice  to  the 
thick  ice  categories.  Consequently,  the  large  amount  of  fall  and  winter  ridging  shown  in  Fig¬ 
ure  24  indicates  that  a  large  amount  of  open  water  is  also  being  created  then.  With  regard  to 
decay,  on  day  180  the  standard  case  has  25%  open  water,  as  contrasted  to  no  open  water  for 
the  thermodynamic  case.  Also,  the  net  annual  export  of  ice  from  the  basin  amounts  to  about 
0.09  Sv,  a  value  close  to  that  estimated  by  Aagaard  and  Greisman  (1975).  For  comparison, 
the  net  outflow  in  the  high-strength  case  amounts  to  only  about  0.045  Sv.  Note,  however, 
that  the  outflow  rates  do  fluctuate  substantially  (Fig.  21),  a  feature  that  is  commensurate 
with  observations  by  Vijne  (1976). 

In  general.  Figure  24  illustrates  the  important  role  ridging  plays  in  the  mass  balance.  On 
the  average,  ridging  transfers  thin  and  intermediate  ice  to  thicker  categories,  thus  making 
room  for  thinner  ice  and,  hence,  greater  growth.  This  transfer,  in  turn,  offsets  the  ablation 
of  the  thicker  ice.  Note  that  there  is  an  annual  net  ablation  of  thick  and  intermediate  ice.  By 
analyzing  the  transfer  by  growth  as  well  as  by  ridging,  it  can  be  shown  that  the  net  intermedi¬ 
ate  ice  melt  is  partially  offset  by  increased  vertical  growth  as  well  as  by  outflow  and  ridging. 
However,  the  thick  ice  net  melt  is  almost  all  balanced  by  ridging  and  outflow. 

While  there  are  certain  quantitative  differences,  these  mass  balance  results  are  qualitative¬ 
ly  similar  to  observational  estimates  by  Koerner  (1973).  Based  mostly  on  ground  observa¬ 
tions  in  the  central  basin,  Koerner  (1973)  deduced  a  total  annual  ice  growth  of  about  1.1m, 
peaking  in  a  thickness  of  about  3.7  m.  He  estimated  this  growth  to  be  balanced  by  about  0.7 
m  of  outflow  and  0.4  m  of  ablation.  Koerner  also  determined  that  about  40%  of  the  ice  grew 
over  ice  less  than  1  m  thick,  and  that  0.20  m  of  the  new  ice  ended  up  in  ridges.  The  simula¬ 
tions  performed  here  yield  a  similar  dominance  of  growth  by  thin  ice  (and  concomitant  in¬ 
corporation  of  new  ice  in  ridges)  but  yield  a  smaller  maximum  ice  thickness.  Also,  there  is 
substantially  more  ablation  and  melting.  Part  of  this  difference  is  likely  from  a  more  com¬ 
plete  inclusion  of  peripheral  areas  of  the  basin  in  the  simulated  estimates.  However,  it  may 
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time-dependent  Maykut-Untersteiner  (1971)  model  and  includes  a  seasonal  snow  cover.  In 
addition,  larger  values  for  the  sensible  and  latent  heat  bulk  aerodynamic  coefficients  are 
used.  (The  smaller  values  used  by  Parkinson  and  Washington  [1979],  and  hence  here,  can  be 
justified  because  geostrophic  rather  than  surface  wind  speeds  are  being  input.) 

In  the  February  case  the  dominant  characteristic  in  the  heat  balance  is  a  rapid  decrease  in 
the  sensible  and  latent  heat  components  as  the  ice  thickens.  The  ice  sensible  heat  losses  simu¬ 
lated  here  are  somewhat  larger  than  Maykut's  values,  primarily  due  to  a  larger  amount  of  in¬ 
coming  longwave  radiation.  (The  difference  in  incoming  longwave  radiation  may  be  deduced 
by  examining  the  open  water  case.)  The  larger  radiation  used  here  tends  to  give  a  higher  sur¬ 
face  ice  temperature  and,  hence,  a  greater  sensible  heat  loss.  This  is  especially  pronounced  in 
the  case  of  thick  ice,  where  Maykut’s  smaller  longwave  radiation  allows  the  surface 
temperature  of  the  snow-covered  ice  to  drop  below  the  air  temperature,  and,  hence,  change 
the  sign  of  the  sensible  heat  transfer. 

In  the  growth-melt  transition  case,  the  main  difference  is  in  the  net  shortwave  radiation 
terms.  The  smaller  thin  ice  albedos  used  by  Maykut  give  a  greater  absorption  of  shortwave 
radiation.  This  in  turn  results  in  a  greater  surface  temperature  and  larger  sensible  heat  flux 
than  in  the  simulation  performed  here.  This  condition  is,  however,  reversed  in  the  very  thick 
ice,  where  a  Maykut  snow  albedo  greater  than  0.7S  allows  less  shortwave  radiation  to  be  ab¬ 
sorbed.  Note  also  that  in  this  transition  period  the  sensible  heat  loss  is  less  dependent  upon 
thickness.  This  arises  from  the  smallness  of  the  conductivity  term  (which  characterizes  the  ice 
thickness),  ascribable  to  a  small  temperature  differential. 

It  is  notable  that  despite  differences  in  parameterization  and  forcings,  the  simulated  ice 
growth  rates  are  relatively  similar  in  both  the  February  and  June  cases.  This  demonstrates 
the  negative  feedback  character  of  the  heat  budget  in  the  presence  of  an  ice  cover.  Basically, 
the  system  attempts  to  minimize  the  effect  of  changes  of  one  component  by  the  adjustment 
of  other  components.  Larger  incoming  longwave  radiation  will,  for  example,  be  offset  par¬ 
tially  by  a  greater  surface  temperature  (which  in  turn  causes  more  outgoing  radiation)  and 
F*  ially  by  an  increased  sensible  heat  loss.  In  the  open  water  case  these  feedbacks  cannot 
operate  since  the  surface  temperature  is  fixed. 

Finally,  in  the  26  June  case  the  surface  temperature  in  all  the  ice  cases  is  fixed  at  freezing. 
As  a  consequence  the  ice  melting  rate  is  independent  of  thickness,  since  all  snow-free  ice  al¬ 
bedos  are  assumed  equal.  There  is  a  difference  in  conduction,  with  the  thinner  ice  having 
slightly  greater  conduction  of  heat  into  the  mixed  layer  and,  hence,  more  bottom  melting. 
Note,  however,  that  in  all  cases  most  of  the  melting  is  surface  ablation.  The  greater  heat  ab¬ 
sorption  of  the  open  water  arises  from  its  lower  albedo. 


APPENDIX  C:  THICKNESS  FINITE  DIFFERENCE  CODE 


Thickness  partition 

For  the  thickness  differencing,  an  arbitrary  number  of  irregularly  spaced  thickness  levels 
are  employed  with  a  center  thickness  A,  and  thickness  limits  hi  to  hi,.  Open  water  is  con¬ 
sidered  simply  by  having  one  thickness  category  centered  at  zero  thickness.  In  order  to  avoid 
discontinuities,  uie  thickness  mesh  is  forced  to  vary  smoothly.  A  flexible  procedure  used  in 
certain  numerical  ocean  models  (vis  Bryan  1969)  is  to  allow  the  spacing  between  categories  to 
vary  according  to  a  Gaussian  distribution.*  Using  this  procedure,  one  obtains  a  partition  of 
center  thicknesses  by 


hm 


0  for  m  =  1 

/»„-i+Ci+c2U-expt-(/n-l)Vc,])  for  m  *  1 


(Cl) 


where  cw  c2  and  c,  are  constants.  Thickness  limits  are  then  determined  by 

I  (ht  -  Ci)/2  for  m  =  1 

(C2) 

Ih^-i-hJ.-i  for  m  >  1. 

Clearly,  by  varying  c,,  c2  and  c,  a  wide  variety  of  different  thickness  partitions  may  be  ob¬ 
tained.  The  particular  arrangement  used  in  these  simulations  (see  Fig.  1)  was  by  using  c,  = 
0.2  m,  c2  =  7.0  m,  and  c,  =  110. 

By  use  of  this  thickness  partition,  the  areal  thickness  distribution  is  defined  at  the  interior 
of  the  thickness  interval  A/to  A,i,  and  is  denoted  by  g,.  For  computation,  an  auxiliary  dimen¬ 
sionless  distribution  function  gi  is  defined  by 


g,  =  dA,  g. 


where  AA,  =  A,i,  -hi.  For  future  differencing,  g,  is  considered  to  have  a  mean  thickness  of 
h,.  The  definition  of  g,  is  that  it  is  the  fraction  of  area  covered  by  ice  with  thickness  between 
A/and  A/,.  The  growth  rates  for  each  category  are  also  defined  at  the  interior  of  the  thickness 
category  and  are  denoted  by/  =  f(hl). 

Thickness  finite  differencing 

To  carry  out  the  operations  in  the  time  marching  equations  (eq  10- 15)  it  is  necessary  to 
have  thickness  finite  differences  for  the  term  (fg)„,  where  the  subscript  A  denotes  a  thickness 
finite  difference.  In  addition,  a  finite  difference  procedure  is  needed  for  the  redistribution 
and  the  lateral  melting.  The  finite  differences  for  the  lateral  melting  steps  are  relatively  trivi¬ 
al.  In  particular,  integrals  over  g  are  simply  replaced  by  sums  over  g,.  However,  the  other 
terms  are  more  complex.  Consequently,  the  differencing  schemes  for  the  (/g)»  term  and  the 
redistribution  are  described  here.  The  other  spatial  finite  differences  are  identical  to  those 
used  by  Hibler  (1979)  and  are  described  there. 


*  Personal  communication  with  M.  Cox,  Oeophysical  Fluid  Dynamics  Laboratory,  Princeton,  New  Jersey,  Novem¬ 
ber  1979. 
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There  are  two  constraints  that  must  be  preserved  by  the  finite  difference  (/#)*.  One  is  that 
the  sum  of  the  changes  in  g  are  zero.  This  follows  from  the  fact  that  vertical  growth  does  not 
change  the  total  area  of  ice  and  open  water.  The  second  is  that  the  new  total  thickness  must 
differ  from  the  old  thickness  by  an  amount  equal  to  the  ice  grown  (or  melted).  These  con¬ 
straints  can  be  expressed  in  finite  difference  form  by  (where  M  thickness  levels  are  assumed) 


I(/DJ/  =  0  (C3) 


M 

l tfftl'A.--  ^  Mo  (C4) 

/=! 

A  finite  differencing  satisfying  these  constraints  may  be  obtained  by  using  a  form  of  up¬ 
stream  differencing  (see,  e.g.,  Mesinger  and  Arakawa  1976).  In  particular  a  flux  (F)  at  the 
thickness  boundary  hi  for  /'  =  2  to  M  is  defined  by 

F  =  [Max(0,  /,.,)  g,.,  +  Min(0 ,/,)  g,]  /(F-F-,). 

To  ensure  conditions  C3  and  C4,  it  is  further  insisted  that 

F  =  =  0. 

Using  these  fluxes,  the  finite  differences  for  the  (/g)»  term  are 

K/kM.  = 

Multiplying  both  sides  by  Ah  =  hi,  -  hi  gives 


M 

E 

/=! 


[(/g)»]<  =  F-i  -F.  (C5) 

The  condition  of  eq  C3  is  satisfied  since  the  FIs  cancel  in  pairs  in  the  sum.  To  see  that  eq  C4 
is  satisfied,  consider  an  arbitrary  growth  rate,  say/,  *0.  Then  the  left-hand  side  of  eq  C4  be¬ 
comes 


~  Fjhj  +  hj-iFj-, 


fjgAhj-i  -  F) 
F  ~  F- 1 


-fjgj- 


The  exception  to  this  occurs  when  open  water  is  melting  or  the  thickest  category  is  freezing. 
In  these  cases  the  conditions  F  =  Fw„  =  0  prevent  terms  necessary  for  satisfying  eq  C4 
from  being  included.  However,  in  the  case  of  open  water  melting  this  heat  will  be  correctly 
accounted  for  by  additional  bottom  or  lateral  melting  as  described  earlier  (see  eq  3).  To  cor¬ 
rect  the  thick  ice  freezing  case,  any  growth  over  thick  ice  is  reapportioned  over  the  other 
growth  rates  in  a  manner  similar  to  eq  3.  In  practice  this  amount  can  be  made  very  small  by 
taking  the  thick  ice  category  to  be  very  thick. 

With  this  particular  finite  difference  code  (and  dynamical  code  described  previously)  the 
essential  stability  requirement  is  a  Courant-Friedrichs-Lewy  condition  for  the  advection 
terms,  namely 
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for  the  spatial  advection  and 


A/s  |<A,-A,-,)j  minffMa x(0, /,-,)]-,  -[Min(0,/)]-'), 

for  the  thickness  advection.  Of  these  two  terms  the  thickness  advection  term  will  normally  be 
most  critical. 

Mechanical  redistribution 

Three  mechanical  redistribution  operations  are  required  at  each  time  step  in  the  numerical 
scheme.  The  operations  consist  of  creating  open  water,  carrying  out  an  ice  thickness  redistri¬ 
bution  and  estimating  the  ice  strength  at  the  end  of  the  time  step.  For  the  open  water  crea¬ 
tion,  all  that  are  needed  are  finite  differences  for  the  strain  rate  tensor.  These  are  calculated 
in  a  way  identical  to  that  used  in  the  dynamical  code  and  are  described  in  Hibler  (1979).  (This 
reference  also  contains  explicit  expressions  for  the  stress  state  atJ  in  terms  of  the  strain  rates 
and  ice  strength.)  For  the  other  two  operations  it  is  necessary  to  introduce  finite  differences 
for  redistribution.  For  this  purpose  it  is  useful  to  define  a  discrete  ridging  operation. 

Using  the  thickness  partition  notation  introduced  in  Appendix  A,  one  defines  a  discrete  re- 
distributor  7  by 


yihfhj)  = 


y{h„h„)  = 


j"  7 (hi,  h')dh'  for  j  <  M 

V 


w 

/  7  (hith')dh'. 


To  ensure  that  7  is  normalized  it  is  insisted  that  M 
M 

T.  y{.hj,h)h(  =  hj. 


In  addition  to  7  it  is  useful  to  define  an  auxiliary  function  X,  by 


\  =  2J  ****'>• 


Essentially,  X,  is  the  total  reduction  of  area  when  a  unit  area  of  ice  of  thickness  h,  is  ridged. 
To  demonstrate  how  a  redistribution  is  carried  out,  consider  that  an  area  reduction  of 

M 

**  -  E  *_1 

/-1 

is  required.  To  redistribute  ice  an  initial  set  of  category  reductions  (denoted  by  g7)  are  esti¬ 
mated  by 


In  this  integral  gt — in  the  definition  of  P(h) — is  considered  to  be  constant  over  the  interval 
h/+  hU |.  Defining  AA '  by 


M 

AA'  =  ^  X,  g? 

i=l 

a  corrected  set  of  category  reductions  are  calculated  by 
f?*+(^L4/A/4 ')  . 

The  new  g,’s  are  then  given  by 
M 

£  Khj.hdgT. 

y-i 

To  estimate  strengths,  the  simplest  approach  is  to  artificially  enlarge  all  g  values  by,  say, 
0.1%.  The  potential  energy  is  then  calculated,  a  redistribution  carried  out,  and  a  new  poten¬ 
tial  energy  calculated.  The  difference  in  potential  energies  divided  by  the  strain  then  yields 
the  ice  strength. 

Formally,  this  procedure  consists  of  first  defining 

g?  =  1.001  g,. 

As  in  Rothrock  (197S),  the  initial  potential  energy  is  given  by 

M 

S?  =  C  ^  g !*  h}. 

/=i 


After  redistribution  g?  becomes  g,  with  potential  energy 
M 

S,  =  C  ^  g,  hi 

trf 


Insisting  that  the  compressive  strength  times  the  divergence  rate  equals  the  rate  of  potential 
energy  change,  one  obtains  the  following  equation  for  P*: 

p.  (0-000  _  (S,-S,«) 

P  At  At 

In  these  expressions  C  is  set  equal  to  2C»  to  allow  for  frictional  losses. 


A  facsimile  catalog  card  in  Library  of  Congress  MARC 
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